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Designed to serve as a guide for integrating 
interactive problem-solving or simulating computers into a 
college-level physics course, this anthology contains nine articles 
each of which includes an introduction, a student manual, and a 
teacner’s guide. Among areas covered in the articles are the 
computerized reduction of data to a Gaussian distribution, 
"free-style" input which frees the student from format restrictions, 
the e^rects of integrating physics concepts with computer programing, 
harmonic motion, relativistic and nonrelati vistic two-body" 
collisions, simulation of relativistic collisions in a bubble 
chamber, further computer simulations involving mass spectrometers 
and physical systems (including the roulette wheel or Monte Carlo 
approach) , vacuum and low— velocity ballistics in a conversational 
anoroacn , and computerized curriculum open— ending to provide 
increased realism and relevance. The student manual and teacher’s 
guide sections clarify and expand the information presented in each 
introduction. Computer listings, subroutine, and structure maos are 
provided also f or many of the articles. (Author/5?) 
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PREFACE 



The articles contained herein are the outgrowth of a compu- 
ter workshop organized by Dr. John W. Robson* for the Commission 
on College Physics in 1968. He succeeded in forming interested 
individuals into a Computer Working Group which met several times 
during the year; their contributions constitute this volume. 

Almost every contributor presented the Editor with a wealth 
of material and those segments were selected for inclusion which 
seemed to the latter's judgment to fit most harmoniously together. 
This was done with considerable anguish and soul-searching; one 
result of this synthesis is that a given contribution does not 
necessarily represent the author's best or most unique work, but, 
in most cases, only the tip of the iceberg. In consequence, the 
Editor has thumpingly urged the authors of the lost masterpieces 
to make them available to the general public through the Computer 
Library for Instruction in Physics (see the American Journal of 
Physics 35 , 273 (1967)). 

In using the computer in physics education, as distinct from 
research, we are not only interested in ways of solving problems, 
but in how the computer can add a new dimension to the nature and 
content of the curriculum through its influence on the topics se- 
lected and their mode of presentation. In order to remain relevant 
to future needs, basic undergraduate courses must be appropriate- 
ly modified to reflect the new points of view associated with com- 
puter applications, numerical analysis must be integrated into 
course work, and students should be given programming instruction 
at an early stage in their education. 

The use to which we put the computer depends on the available 
facilities . From a pedagogical standpoint these are of three 
types — interactive, semi-interactive and noninteractive. Large, 
expensive installations which process programs in batches are non- 
interactive. Smaller and cheaper machines may allow a student to 
receive his output within a few minutes, alter his program and 
resubmit it, if necessary. These are semi-interactive in that re- 
ceipt of output is immediate, but errors or program alterations 
necessitate terminating the program. The future widespread use of 
the computer in the physics curriculum will result from the avail- 
ability of interactive terminals, where the student can receive 
output immediately, correct errors and input new instructions or 
data without terminating his connection with the computer. Hence, 
the emphasis in this volume has been on interactive and semi-in- 
teractive facilities. 

The computer may be used in the classroom in any of four ways: 



*Presently, Dr. Robson is with the Department of Physics at 
the University of Arizona, Tucson. 
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as an administrator, tutor, simulator or calculator. In the first 
mode it may be used simply to administer and grade exams and per- 
form other onerous clerical chores . In the second mode it may , in 
addition, tutor the student; i.e., correct his errors by hint, 
precept or example and so lead him along various paths which are 
determined by his previous answers and designed to aid him to over- 
come his deficiencies . 

The computer can simulate physical reality either as a "black 
box" or as a loaded roulette wheel. The "black box" may be a pro- 
gram representing a physical system into which the student enters 
the values of certain physical quantities and then observes the 
output from the program in place of an actual experiment. On a 
time-shared teletype the student can easily "fiddle with the knobs," 
i .e . , change his inputs and observe the results . In the roulette 
wheel, or Monte Carlo, method the computer generates random numbers 
which can be used to simulate phenomena in which chance in a factor . 
For example, one could simulate a baseball game using the batting 
averages of the players. A player batting .325 would be allowed a 
hit each time the three-digit random number representing his time 
at bat was between 000 and 324 and called "out" if it was between 
325 and 999. Such methods can be applied to problems of gaseous 
diffusion, radioactive decay, scattering, etc. 

However, it is as a calculator and solver of problems that the 
computer should have its greatest impact on physics education. And 
rather than have students use programs they do not understand, it 
would be preferable to integrate the computer into physics at an 
early stage. The major problem is the lack of textual materials 
and programs, and the need for wider dissemination of those which 
do exist. It is to serve this need that this volume was conceived, 
and for this reason, also, programs are presented in FORTRAN or 
BASIC, the two most popular languages in use today. 

Since the individual authors generally performed their work 
prior to the organization of the Computer Working Group, some dup- 
lication was inevitable, and, in the case of the harmonic oscilla- 
tor — probably the world’s most "programmable" elementary problem — 
was considered beneficial. Thus, we have presented the work of 
Vierling which primarily illustrates the application of an advanced 
fourth-order Runge-Kutta method to the theoretical problem, as well 
as that of Grimsrud, designed for use in connection with an ele- 
mentary pendulum experiment. 

The volume opens with "Data Reduction" by Smith, which is in 
the nature of a prerequisite to computer-oriented physics, since 
it deals with the use of the computer to reduce data to a Gaussian 
distribution and also describes an auxiliary program which may be 
used to interpret "free— style" input, in which the students* input 
to the program is freed from the usual tedious and confusing re- 
strictions on format. Following Vierling and Grimsrud is Winder’s 
article, which describes relativistic two-body collisions, Grimsrud 
having already presented the nonrelativistic case. This leads, 
quite naturally, to Mikkelson’s simulation of relativistic colli - 
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sions in a bubble chamber. The next articles, by Harbron & Miller 
and Kinsey & Kenyon, further illustrate the use of the computer 
as a simulator; the latter illustrating, too, the Monte Carlo ap- 
proach . 

Although we have generally elected to avoid the tutorial mode 
in this monograph, Jalbert's treatment of vacuum and low-velocity 
ballistics does provide an example of a conversational approach. 
Finally, Blum's article attempts to demonstrate how the computer 
can be employed in the construction of an open-ended curriculum in 
the sense that without much broadening of the conceptual base, as 
given by Jalbert, the student can be enabled to treat problems of 
greatly increased realism and relevance. 

Although the Editor has required contributing authors to cast 
their articles in the same format — Introduction, Student Manual, 
and Teacher's Guide — a conscious attempt has been made to preserve 
individual nuances of style and approach on the grounds that these 
subtleties are themselves of interest to the practicing pedagogue. 
However, the Editor takes no responsibility for the absolute rec- 
titude of the contributing authors, wich the exception of the fi- 
nal article in this anthology. Criticisms and suggestions are 
welcomed, and remarks directed to particular authors will be for- 
warded to them by the Editor. Individual articles may be ordered 
separately through the Commission on College Physics, hence the 
unusually verbose footnotes scattered through the text. 

The Editor acknowledges the invaluable assistance of Mrs . 

Faye von Limbach who prepared the typescript and the flow charts 
and whose pithy observations often served as a useful stimulus . 
Equally valuable were the tireless and enthusiastic efforts of Mr. 
Lee A. Fowler who checked out many of the programs in this work. 
Miss Kathryn E. Mervine also assisted in the editing, and the 
project enjoyed the wholehearted support of Dr. John M. Fowler, 
Director of the Commission on College Physics . 

Members and friends of the Computer Working Group, in addi- 
tion to those mentioned above, included Alfred Bork, University of 
California, Irvine; David J. Cowan and Richard T. Mara, Gettysburg 
College; S.A. Elder, U.S. Naval Academy; Russell K. Hobbie, Uni- 
versity of Minnesota; Arthur Luehrmann, Dartmouth College; Anatole 
Shapiro, Brown University; Harold Weinstock, Illinois Institute of 
Technology; Ronald Winters, Denison University; Claude Kacser, 
Leonard Rodberg and Sanders N. Wall, University of Maryland. 

By and large the papers presented here do not claim to be 
unique or particularly efficient. Hopefully, as a collection, they 
are all of pedagogical interest as regards their presentation, the 
problem treated and the experiences of their creators . We hope it 
may serve as a guide and companion to those who wish to integrate 
computers meaningfully into the context of their course work. 

Ronald Blum 

Commission on College Physics 
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INTRODUCTION 



This paper describes "GAUS," a program designed to compute 
the mean and standard deviation of a collection of laboratory data, 
as well as its goodness of fit to a Gaussian (normal) distribution 
by means of the chi-squared test. Also described is a general 
purpose auxiliary program, "FREE," which allows the students to 
input data in a format- free mode. This greatly simplifies the use 
of the computer in the classroom or laboratory, since data can be 
written in a natural way and the students need not learn the in- 
tricate and often tedious details of formatting, a fruitful source 
of time-consuming errors. 

Use of computer facilities to support undergraduate physics 
laboratory instruction started at Coe College in the academic year 
1965-66. A number of programs were written to facilitate data re- 
duction in an intermediate electricity and magnetism laboratory 
and later "GAUS" was developed to help students in an introductory 
laboratory gain insight into error of measurement. For the first 
two years Coe made use of the IBM 7044 located at the University 
of Iowa, twenty miles distant, via a courier service that provided 
overnight turn-around. Initially, "GAUS" was used via optical 
sense-mark cards developed by the Measurement Research Center in 
Iowa City and requiring only a number two pencil for marking. 

In the spring of 1967, under an NSF grant, Coe installed an 
IBM 1130 with 8K core and 1/2 M disk; that summer existing programs 
were modified for the 1130 and, with the help of high school stu- 
dents, an extensive library of simple but generally useful subrou- 
tines was developed to simplify further programming. This included 
the first version of FREE and its associated supportive subroutines. 
In the fall of 1967 and again in 1968, the author taught an intro- 
ductory physics laboratory with calculus, in which computer use 
played a major role. The programs were introduced to the students 
as black boxes with little or no reference to their details of op- 
eration. Emphasis was on least squares fits of experimental data 
for linear acceleration, damped linear acceleration, and damped 
simple harmonic motion. It has taken typical students three to 
five weeks to really grasp the concept of least squares fits in 
these situations. About 280 students have made use of "GAUS" and 
found it a stimulating educational experience. 

On the basis of three years' experience, the author has the 
following recommendations to make: 

1. Simplify input-output problems with a free style card 
reader. 

2. Develop a balanced diet of very simple programs that 
students can write and/or modify, along with complica- 
ted programs such as those used here. 
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3. Dor? • t expect to teach students very much FORTRAN while 
they are simultaneously engaged in more traditional 
laboratory experiments. 

4. Be alert for students who have "personal relations" 
problems with the computer. Some may not get assigned 
work done because they feel an antipathy toward the 
computer, others may lag behind due to their infatua- 
tion with it and associated diversions. Both will need 
sympathetic help! 

The author wishes to acknowledge the assistance of Dr. Joseph 
Kasper of Coe College and Dr. Ronald Blum of the Commission on 
College Physics, who are responsible for much of the documentation 
on the Gaussian distribution in the Student Manual which follows. 
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STUDENT MANUAL 



Gaussian Distribution 

Suppose that a measurement of a physical quantity is made, 
where it is known that completely random deviations in the mea- 
surement occur. Suppose further that the measurement is made 
repeatedly, with all sources of systematic error eliminated. 

The theory of probability demonstrates that if one plots the 
frequency of occurrence of a given value of N as a function of 
the value of N, and if one collects very large numbers of data, 
then the graphical representation of the results will ideally be 
in the form of a symmetrical curve known as a Gaussian curve, or 
Gaussian distribution. 

A typical Gaussian distribution, symmetric about the mean 
value, N, is shown in Figure 1; the ordinate, P(N), is the rela- 
tive probability (or "probability density") that any particular 
measurement, N, will occur. While in reality the quantity mea- 
sured may only take on certain discrete values, we shall find 
that the continuous distribution is in fact a very useful tool 
for the reduction of such data. The analytical formula for this 
curve is __ 

P (N) = (2tt a 2 )“ V2 exp [- (N - N) 2 /2 a 2 ] (1) 

where a 2 is an independent parameter known as the variance of 
N and its root, cr, is called the standard deviation of N; exp (x) 
= e x , another notation for the exponential function. The constant 
multiplier (2ira 2 )“^2 is chosen such that the total area under the 
curve between the limits - » < N < 00 is exactly unity. Thus, it 
must follow that the area under the relative probability curve 
between, say, N^ and N 2 is^ the theoretical probability that a 
measurement of N will fall between those two values. Further- 
more, although there is some inconsistency in special cases where 
N can never actually be negative, such as the count rate from a 
radioactive source, this discrepancy is negligible the area 
under P(N) for negative N is only of the order of a few per- 
cent, which is generally the case in our experiments. 

_ Figure 1 shows two Gaussian curves with the same mean value 
N = 9, but two different values of standard deviation, o. Al- 
though both curves peak at N = N = 9, the curve with greater de- 
viation is wider and, consequently, lower, since the area under 
each curve must be unity. If one made but one measurement of the 
physical quantity involved, and knew that the lower curve applied 
then he would have little confidence in it , because on the grounds 
of probability there would be a good chance that his one value was 
considerably greater or smaller than the mean. However, if he 
made one measurement and knew that the higher curve applied, his 
confidence in the measurement would be much higher. This degree 
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of confidence, or probable deviation from the mean, is stated in 
various ways . 




Figure 1. 



Gaussian distributions for N = 9, a 



1, and 0 = 3. 



It is a truly remarkable aspect of physical measurements in- 
volving random errors that: 

a) If a measurement of some physical quantity is made, 
then, because of random fluctuations, it cannot be 
taken as entirely reliable, even if there are no 
instrumental or human errors present. However, the 
known Gaussian distribution can be used to furnish 
a meaningful estimate of the validity of the mea- 
surement. 

b) If sets of measurements are made, then they can be 
compared with the ideal Gaussian curve, and devia- 
tions in the actual distribution can be used to 
check the validity of the results . 

The standard deviation is a measure of the random variability of 
individual measurements and hence of the confidence we may place in 
them. It is often useful to think in terms of units of o about 
the mean value, as shown in Figure 2. One can see that 68.3% of 
the time an individual measurement taken at random will lie within 
± 0 of the mean value, and 95.7% of the time within ± 20 of the 

mean. Furthermore, in statistical theory, it is shown that the 
mean value, while it may not be precisely the same as the true val- 
ue, is the best estimate we can form of the true value; our results 
are no less "scientific" for having taken account of the harsh re- 
alities of random errors. It should also be clear from Figure 2 
that when N is greater than 2 o the fact that P(N) is also de- 
fined for N < 0 will not cause any serious errors in estimating 
probabilities in radioactive counting experiments from the Gaussian 
distribution . 

In general N and o are independent of each other; while 
repeated measurements of the period of a pendulum yield an average 
value depending primarily on the physics of the pendulum, the stan- 



6 



dard deviation of the measurements depends on the instrument and 
technique of measurement. However, there is one important situa- 
tion in which a very interesting relationship exists between N 
and 0 . Whenever the measurement consists of counting independent 
events occurring at random in successive intervals of time (or 
space) , then 



provided that N is moderately large (i.e., N > 10). 

For example, if one counted events with a Geiger tube over 
a radioactive source and got 100 counts in some time interval, 
such as ten seconds — hence, the standard deviation is o = 10 
counts — then 68.3% of the time such a single measurement would 
be within ± 10 counts of the average number of events due to 
that particular source. That is, were we to estimate that 90 < 
N < 110, we would be right 68.3% of the time, or in 683 cases 
out of 1,000 such measurements. The counting rate, R, would be 



with a "confidence level" of 68.3%. Were we willing to make a 
looser prediction we could say R = (100 ± 20) counts/ (10 sec), 
with a confidence level of 95.7%. 



The standard deviation also has the significance that if we 
form the quantity (N - N) 2 for each measurement of N and average 
all such values, we obtain a 2 ; hence, a may also be known as the 
root-mean-square (rms) error of a variable, in counting-situations. 



estimating 0 . However, in most situations one needs several mea- 
surements to form an accurate estimate of the standard deviation, 
according to the formula 



0 






R = (100 ± 10) counts/ (10 sec) 
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Figure 2. Gaussian distributions in intervals of 0 



where 0 = /jj , even one measurement of N affords a basis for 
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where a total of M measurements of counts have been made; is 

the result of the i-th measurement. (It is actually most correct 
to use M-l in the denominator of the above expression; however, 
this is a subtle theoretical point, and if M is large this is a 
negligible effect.) 

Now suppose we have collected some statistical data in the 
form of repeated measurements . What can be made of them? They 
can be compared with the expected Gaussian distribution 
by making what is called a "chi-squared test" . The quantity "chi- 
squared" is computed by first preparing a histogram, or bar-graph, 
showing the number of times ( frequency ) each measurement appeared, 
as a function of the value of that measurement. Actually, each 
histogram bar represents the number of times that measurements 
lead to values falling in a certain interval;, e.g., one bar might 
.represent how many times the count was between 90 and 100, the 
next bar how many times the count was between 100 and 110, etc. 

When the completed histogram is at hand, it can be superim- 
posed on the ideal or expected Gau.ssian distribution. There will 
be a difference between the actual number of values, nj, observed 
to fall in the j-th interval and the e xpected number of values, mj, 
which the Gaussian distribution predicts will fall in the interval. 
If, for each such interval, we compute the quantity (nj-mj) 2 /mj, 
the sum over all J intervals is called the chi-squared statistic. 



X 
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j V J ( n J ~ mj) 2 
j=i m- } 



The better the fit of distribution to data, the smaller is x 2 - 



Probability theory tells us, for randomly distributed errors, 
the probabilities of obtaining different values of x 2 f or differ- 
ent values of J;* thus, one can check equipment or data for sys- 
tematic or nonrandom errors by comparing the fit of the data to a 
Gaussian curve. This information is commonly found in handbooks 
in the form of tables which give the probability that x 2 will 
equal or exceed a certain numerical value. A sample of one line 
from such a table is shown below: 



Number of Intervals, J There is a probability of 

0.99 0.90 0.50 0.10 0.01 

. that the calculated x 2 is equal to or greater than 



19 



7.63 11.65 18.34 27.20 36.19 



Thus, there is no unique answer to the question: When is a 

fit good or bad? Instead, there is only a probabilistic answer. 
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For example, as the excerpt from the table shows, if one measured 
a count rate 19 times, plotted the distribution, plotted the ex- 
pected Gaussian distribution, computed the value of x 2 / and found 
X 2 to be 10, the fit would be quite good. Specifically, there 
would be more than a 90% chance that if the experiment were repeat- 
ed the x 2 would be greater than 10. On the other hand, if x 2 
turned out to be as high as 40, the fit would be considered very 
bad and it would be unlikely that the distribution represents a 
random selection from a set of Gaussian distributions. 

In collecting data by measuring some quantity over and over, 
peculiar instances sometimes arise. For example, if one were con- 
sistently getting between 360 and 440 counts in one-minute inter- 
vals from a radioactive source, and then in a given one-minute in- 
terval got only 120 counts, one would be suspicious of that result. 
According to the Gaussian curve, such an anomalous count is per- 
fectly possible, for the Gaussian curve runs to infinity in both 
directions from the mean. However, the probability of counts far 
from the mean drops off more and more rapidly, the farther they 
are from the mean. If one gets a count of 120 in a one-minute 
interval, after consistently getting close to 400 each minute, the 
thing that should really disturb us is not that this is totally 
impossible, but that it is highly unlikely . In fact, it may seem 
so unlikely that one would not want to include the result at all, 
because it is not typical and would throw the mean value off with 
more weight than it deserves. In short, one is tempted to reject 
that far-off value. 

However, such subjective selectivity constitutes rather wain- 
ton tampering with the scientific data. The prime requisite for 
scientific honesty and objectivity is to let nature speak for her- 
self, rather than to interpose the subjective bias of the experi- 
menter. The solution to this problem is to adopt a specific cri- 
terion, expressly stated, for the acceptability of data. One such 
is "Chauvenet ' s Criterion" which states that an observation should 
be discarded if the probability of its occurrence in the set of 
observations is equal to or less than 1/ (2K) , where K is the 
number of observations. The table below, assuming a Gaussian dis- 
tribution, gives the maximum acceptable departure of any one 
reading from the mean in units of a : 

Number of Observations Maximum Departure Acceptable in Units of Standard Deviation 
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1.65 


10 


1.96 


20 


2.24 


50 


2.58 


100 


2.81 


200 


3.02 


500 


3.29 


1000 


3.48 



Interested students may find it instructive to investigate 
the Poisson distribution, which, when its mean is large, has the 
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shape of the Gaussian distribution. The ideal Poisson distribution 
has the peculiar property 'that, regardless of the magnitude of the 
mean, the standard deviation is equ^.l to the square root of the 
mean of the distribution of values. 



Instructions for Using the Computer Program " GAUS " 

* 

(1) Make repeated measurements of some physical quantity according 
to whatever instructions you have been given for the particular 
laboratory session. Take care to record the values obtained in the 
order you obtain them . 

(2) Punch into standard IBM data cards the values you obtained: 

a. Separate different values by one or more blank positions 
on the card. 

b. Punch decimal points only for values with a decimal 
fraction part. 

c. Punch as many values on one card as you can, then con- 
tinue on another card. 

d. Do not split a given value between two cards. 

e. Punch errors may be "erased" by overpunching each column 
of the given value xvith an X, _if column one is blank . 

f. After your last value, enter "9999 w in your data card, 
no decimal point! 

(3) Make up a data deck as follows using pre-pun ched cards if they 
are supplied; otherwise, punch your own additional cards as needed. 
The symbol "b" does not stand for the letter "B" punched on a data 
card; it stands for a blank column with no hole punched in it. 

Blank columns so designated must be present. The initial "//" must 
be in columns one and two of the cards. 

Here is how your deck of cards should look (top line represents top 
card) : 

//bJOBbT 
// bXEQbGAUS 

A card with your name or identifying nui.iber punched. 
Your data cards — in order — as prepared under (2); 

« be sure you enter 9999 after your last value. 

/ /b*b JOBbEND 

(4) Submit your program; or run it yourself, if so instructed. 



Output from GAUS 

The first printed line is your identification card. 

Values rejected according to the Chauvenet Criterion are noted 
next; one or more values may be rejected at a time. After each 
such occasion of rejection the remaining list is again checked to 
see if any values in it should be rejected according to Chauvenet 1 s 
Criterion. Any value that is exactly zero is automatically reject- 
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ed, as it is likely to be due to a keypunch or card reading error. 
You may enter values arbitrarily close to zero, however; e.g. f 
0.0001 is permitted. 

Next, the values remaining are printed out in the original 
order on the left, along with three additional columns: the "Run- 

ning Average," "Running Sigma," and "Running Sigma of Ave." The 
"Running Average" is the average of the value to its left plus all 
preceding values in the first column; the "Running Sigma" is' the 
standard deviation of those values; "Running Sigma of Ave.” is an 
estimate of how confident you can be in the value of the "Running 
Average" value just to its left. It is, in effect, a prediction 
of the standard deviation you might expect to get upon making a 
list of average values obtained in a manner exactly as the one to 
its left was obtained. That is, if the average value, x, of a 
collection of data, x, ( statistics ) is itself considered as a new 
statistic, the theory of probability predicts that a collection of 
average values derived from data samples of M different measure- 
ments will have a standard deviation <j x related to the deviation 

of a single sne, asurement by 

c*x = o x /M 

The last values in these "Running" lists apply to the full list 
given, except for those values rejected. 

The full list is next divided up into successive subsets to 
illustrate the way in which the average and standard deviation of 
successive subsets fluctuates about the average and standard de- 
viation applicable to the full list. Compare these values to the 
last entries under "Running Average" and "Running Sigma". 

The list of subset averages is now considered as a new list 
of values for which the average and standard deviation is computed. 
The same is done for the list of the standard deviations of the 
subsets, which may, in turn, be considered a statistic and assigned 
an average, cx, and a standard deviation, <5o% . 

The first value under "Their Sigma, " which represents the 
standard deviation of the list of averages of subsets, should be 
compared with the M-th entry under "Running Sigma of Ave." above, 
where M = the number of values in each subset of values . These 
two values should be within 10 to 20 percent of each other if your 
list of values has a true Gaussian shape. 

Your list of values is again printed out, but this time or- 
dered according to value; also listed is the departure of each 
value from the mean value, i.e., the last value in the list "Run- 
ning Average." 

Next, a histogram is printed out along with an indication of 
the interval of values included in each histogram bar. This is 
followed by a numerical comparison of the observed histogram with 
the expected. "CHI2" heads the list of contributions to the chi- 
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squared sum. "Lower" and "Upper" are bounds on the histogram bars 
in terms of departures from the mean. If the data is distributed 
according to a Gaussian distribution, then theory tells us that 
the value of x 2 for data selected at random will be less than 
the number of histogram bars at least 50% of the time, the exact 
figure depending upon the exact number of bars. Thus, as a rule- 
of-thumb, if x 2 is significantly larger than the number of his- 
togram bars, we would have reason to question the validity of the 
Gaussian distribution as applied to that particular set of data. 
(Incidentally, regrouping the data to make more or less bars also 
affects the value of x 2 accordingly — you can't beat the system!) 




TEACHER'S GUIDE 



Approximately 280 students have used GAUS over the past three 
and a half years . These have been students of high school age in 
summer NSF science programs as well as students in our introductory 
physics courses. The student response has been generally quite 
good, with a significant number of students showing genuine enthu- 
siasm. 

A variety of measurement procedures have been used to generate 
data for GAUS : 

1. Repeated measurements of the period of a simple tor- 
sional pendulum made by hanging a rod on magnetic 
tape. 

2 . Repeated measurements of transit time for a car on 
an air track. 

3 . Radioactive decay counts — each student instructed to 
set his counting equipment so as to get some assigned 
average count rate, at least approximately, so the 
class as a whole can check the prediction that stan- 
dard deviations will go as the square root of the 
mean values. 

Class and/or laboratory discussion is directed to the useful- 
ness of the standard deviation associated with measurement. In 
particular this is related to the question of whether certain func- 
tional relations describe the relationship between experimental 
values to within a reasonable degree of accuracy. The utility of 
the Gaussian distribution needs to be emphasized and illustrated 
by actually making use of it in laboratory work. 

In selecting an experimental procedure to generate data for 
GAUS, the following considerations are important: 

1. It should be possible to obtain about 100 indepen- 
dent measurements /hour; GAUS can handle 200 readings 
in its present form. 

2. Repeated measurements should yield a significant 
scatter of values over a continuous range — specific 
values should rarely repeat. The histogram gener- 
ated may look very strange if this condition is not 
met. 

3. Except as you desire it, there should be no trend 
for the repeated measurements to lead to values 
tending to steadily increase or decrease. Note: 
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extended use of a stop watch can lead to fatigue 
with increased response time and variance — this 
might be interesting to study. 

4. Avoid a change in "observer" in the middle of a 
set of measurements , except as you wish students 
to study possible effects of such changes, if any. 

It has been possible at Coe College to have students go dir- 
ectly from the laboratory to keypunches and then to the IBM 1130 
where they could run their data " open -shop " . While this does help 
develop interest, there was also good student response even when 
we had delays due to courier service to and from Iowa City. With 
the IBM 1130 we store the program in "Core Image Format, " so that 
execution starts within about five seconds and runs one to three 
minutes for typical sets of data used. Students are given the op- 
tion of having their decks run for them closed-shop style if time 
or inclination rules against the "hands-on" operation. 

The following pages contain what it is hoped are self-explan- 
atory listings of GAUS and its associated subroutines, followed 
by an Appendix explaining the use of the subroutine for free-style 
input . 
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// * MOD IF 6 8 05 

// DUP 

♦DELETE GAUS 

// FOR 

* IOCS (CARD, 1132 PRINTER) 

♦ONE WORD INTEGERS 

C GAUS USES SUBROUTINES SUBSE, ORDE, COUNT, GEX, CHI, AND FREE 

C GAUS COORDINATES THE USE OF SUBROUTINES ORDE, COUNT, GEX, AND 

C CHI TO MAKE AN ANALYSIS OF AN UNORDERED SET OF VALUES HAVING, 

C PRESUMABLY, A GAUSSIAN DISTRIBUTION. IT ACCEPTS THE K VALUES 

C IN THE VECTOR X(I) AND GENERATES THE FOLLOWING 

C AVE = THE AVERAGE VALUE OF SET OF READINGS 

C S = THE STANDARD DEVIATION OF THE SET OF READINGS 

C INTV= THE NUMBER OF INTERVALS SET UP TO MAKE HISTOGRAM BARS 

C QLX(I) = THE NUMBER OF VALUES IN VARIOUS HISTOGRAM BARS. 

C X(I) = THE READ IN VALUES IN GAUS, BUT COUNT FREQUENCY IN CHI 

C EX (I) = THE NUMBER OF VALUES EXPECTED IN THESE HISTOGRAM BARS. 

C DX(I) = THE DIFFERENCE BETWEEN ACTUAL AND EXPECTED, 2 DIF. USES 

C DX2 (I)= THE SQURE OF DX(I) 

C CHI2 (I) = CONTRIBUTIONS TO THE CHI SQUARED CRITERION PARAMETER 

C QLB (I)= THE LOWER BOUNDS OF THE HISTOGRAM BAR INTERVALS 

C QUB (I)= THE UPPER BOUNDS TO THE HISTOGRAM BAR INTERVALS 

C SX = THE SUM OF THE NUMBER OF VALUES IN THE HISTOGRAM BARS 

C SEX= THE SUM OF THE EXPECTED NUMBER OF VALUES IN THE VARIOUS BARS 

C SDX= THE SUM OF THE ABSOLUTE VALUES OF THE VALUES DX(I) 

C SDX2 = THE SUM OF THE VALUES DX2(I) 

C SCHI2= THE SUM OF THE VALUES CHI2 (I) , THE CHI SQUARED SUM 

C 

C BY PAUL A. SMITH, COE COLLEGE, CEDAR RAPIDS, IOWA 

C PLEASE COMMUNICATE ANY DIFFICULTIES TO THE AUTHOR DIRECTLY 

C 

DIMENSION X (200 ) ,DX(200) ,QLX(18) ,EX(18) ,DX2 (18) ,CHI2 (18) , QLB (18) , 

. QUB (17) , LX (17) , TITLE (40) 

1 FORMAT ( 4 0A2) 

2 FORMAT (1H1, 4 0A2) 

5 READ (2,1) (TITLE (I) ,1=1,40) 

WRITE (3,2) (TITLE (I), 1=1, 40) 

K = 200 
RJT =0.0 
CALL FREE (X,K) 

C TEST TO BE SURE A SET OF VALUES WAS ACTUALLY READ IN 

IF (K) 15,5,15 
15 QK = K 
SX = 0.0 
DO 20 I = 1, K 
20 SX = SX + X(I) 

AVE = SX/QK 
SDX2 =0.0 
DO 40 I = 1, K 
DX (I) = X(I) - AVE 
40 SDX2 = SDX2 + DX(I)**2 
S = SQRT (SDX2/QK) 

C USE LEAST SQUARES FIT FUNCTION REPRESENTATION OF CHAUVENET CRITERION 

C FUNCTION WHICH IS DEPENDENT ON THE NUMBER OF VALUES INVOLVED 

CHAUV = S* (SQRT (8. 39+8. 62*ALOG (QK) /2 . 30259) -2. 08-0. lll*ALOG (QK) / 
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12.30259) 

J = 0 
DO 50 I = If K 

REJECT ALL VALUES EXACTLY EQUAL TO ZERO AS THEY MAY BE KEYPUNCH OR READ 
ERRORS AND WOULD NOT LIKELY BE REJECTED BY CONSIDERATIONS FOLLOWING 
VALUES VERY CLOSE TO ZERO ARE PERMITTED , FOR EXAMPLE 0.0001 
NEGATIVE VALUES ARE PERMITTED 
IF (ABS (X (!' ) -0.0000001) 42,41,41 

REJECT VALUES VERY FAR OUT ON TAILS OF DISTRIBUTION 

41 IF (CHAUV-ABS (DX (I) ) ) 42,42,48 

42 IF(RJT) 43,43,45 

43 WRITE (3/ 44) 

44 FORMAT ( ' 0 SORRY, BUT THESE VALUES LOOK OUT OF PLACE AND SO ARE REJE 
.CTEDV//HX, 'VALUE' 8X, 'AVERAGE' 6 X, 'DEVIATION PERMITTED DEV' ) 

RJT =1.0 

45 WRITE (3, 46) X(I) ,AVE,DX(I) ,CHAUV 

46 FORMAT (' '4F15.4) 

GO TO 50 

48 J = J + 1 
DX(J) = DX (I) 

50 CONTINUE 

END OF EXTREME VALUE REJECTION LOOP 
TEST TO SEE IF ANY VALUES REJECTED 
IF(K-J) 52,60,52 
52 K = J 

IF VALUES WERE REJECTED RECOMPUTE THE RAW SET OF VALUES 
DO 55 I = 1,K 
55 X(I) = DX (I) + AVE 

SINGLE SPACE THE PRINTER 
WRITE (3 46) 

IF VALUES WERE REJECTED RE-CHECK REMAINING LIST FOR EXTREME VALUES 



GO TO 15 

THE RAW SET OF VALUES HAS BEEN CLEANED AND WE HAVE AVERAGE AND STANDARD 
DEVIATION 

MAKE A STUDY OF RUNNING AVERAGE AND STANDARD DEVIATIONS AND 

THEN FOR SUBSETS OF THE FULL SET COMPUTE AVERAGE AND STANDARD DEVIATIONS 

60 CALL SUBSE (X,K) ____ 

PUT THE VALUES AND THEIR DEPARTURE FROM THE AVERAGE IN NUMERICAL ORDER 

CALL ORDE (X, DX, K) 

WRITE (3,62) (X(I), DX(I), 1=1, J) 

62 FORMAT (///'0 ORDERED VALUES DEV FROM AVE '//( ' '2F15.4)) 

COUNT UP HOW MANY VALUES LANDED IN EACH OF CERTAIN HISTOGRAM BAR INTERVALS 
CALL COUNT (K,DX,SP, INTV, LX, QLB,QUB, AVE) 

DO 70 I = 1 , INTV 

70 COMPUTE THE ^EXPECTED GAUSSIAN DISTRIBUTION OF COUNTS FOR HISTOGRAM BARS 

CALL GEX ( INTV , K, EX, SP , S) « TmmT Bn»T«. 

COMPUTE THE CHI SQUARE VALUE FOR OBSERVED VERSUS EXPECTED DISTRIBUTION 

CALL CHI (INTV, QLX, EX, DX, DX2 ,CHI2 , QLB, QUB, SX, SEX,SDX, SDX2 # SCHI2 , 

1AVE , S) 

PERMIT MULTIPLE SETS OF DATA 



GO TO 5 

C WILL CALL EXIT ON ENCOUNTERING A // * MONITOR COMMENT CARD 

END 

// DUP 

* STORE WS UA GAUS 
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SUBROUTINE SUBSE(X,K) 

DIMENSION X(l) 

MAKE STUDY OF RUNNING AVERAGE AND STANDARD DEVIATION OF SET OF VALUES 
DO NOT MAKE STUDY ON A SET OF LESS THAN 9 VALUES 
IF (K-9) 5,10,10 
5 RETURN 

10 WRITE (3,20) X(1),X(1) 

20 FORMATCO OUTPUT FROM SUBROUTINE SUBSET'/' 0 *9X' VALUES ', 3 (8X* RUNNING 
.'//gXf’AS R^AD Iff* SX, 'AVERAGE* 10X, 'SIGMA* ,3X, 'SIGMA OF AVE'/'0'2Fl 
.5.4) 

SUM = X(l) 

Sum 2 = x(l) **2 

COMPUTE RUNNING AVERAGE AND STANDARD DEVIATION FOR SET OF VALUES 
DO 100 I = 2,K 
SUM = SUM + X(I) 

SUM2 = SUM2 + X (I) **2 

AVE = AVERAGE OF SET OF VALUES INCLUDED UP TO THIS POINT 
SDS = STANDARD DEVIATION OF THE SET 
SDA - STANDARD DEVIATION OF THE AVERAGE 
AVE = SUM/FLOAT (I) 

SDS = SQRT { (SUM2-FLOAT (I) *AVE**2)/FLOAT (1-1) ) 

SDA = SQRT ( ( SUM2— FLOAT ( I ) *AVE* * 2 ) /FLOAT ( I* ( I- 1 ) ) ) 

WRITE (3,30) X (I) , AVE, SDS, SDA 
30 FORMAT (’ ’4F15.4) 

100 CONTINUE 
WRITE (3, 101) K 

101 FORMAT ( ' 0 THERE WERE *13,' VALUES IN THE FULL CLEAN SET') 

PREPARE TO MAKE A STUDY OF SUBSETS OF THE FULL SET OF VALUES 

NOT ALL VALUES CAN BE USED IN STUDY OF SUBSETS, SEEK TO WASTE AS FEW 
AS IS POSSIBLE 

COMPUTE APPROPRIATE SIZE OF SUBSETS TO WASTE LEAST NUMBER OF VALUES 
NGK = NUMBER OF GROUPS KEEP 

KQ = NUMBER OF VALUES WASTED IN BEST CHOICE OF NG TO PATE 

N = NUMBER OF VALUES TO BE GROUPED IS EQUAL TO NUMBER IN LIST 

NS = SQUARE ROOT OF NUMBER OF VALUES, ROUGH ESTIMATE OF NUMBER GROUPS 

NI = NUMBER OF VALUES EACH SIDE OF NS TO BE CONSIDERED 

NX = MAXIMUM VALUE TO BE CONSIDERED IN CONSIDERING NUMBER OF GROUPS 

NG = THE NUMBER OF GROUPS CURRENTLY BEING CONTEMPLATED 

NGK = 0 

KG = 1000 

N = K 

NS = SQRT (FLOAT (N) ) 

NI = (FLOAT (N) ) **0.25 
NX = NS + NI 
DO 105 NG = NS, NX 

CHECK TO SEE HOW MANY VALUES WOULD BE WASTED WITH CURRENT NG 
IF (IABS (N- (N/NG) *NG) -KQ) 104,104, 105 

STORE CURRENT BETTER THAN ANY FORMER VALUE OF NG WITH ASSOCIATED KQ 

104 NGK = NG 

KQ = IABS (N- (N/NG) *NG) 

105 CONTINUE 

NR = NUMBER OF READINGS PER GROUP 

NU = NUMBER OF VALUES USED 

SUM = SUM OF INDIVIDUAL VALUES IN SUBSET 

SUM2 = SUM OF SQUARES OF INDIVIDUAL VALUES IN SUBSET 

GSUM = GRAND SUM OF SUBSET AVERAGES 
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GSUM2 = GRAND SUM OF SUBSET AVERAGES SQUARED 

ZSUM = SUM OF STANDARD DEVIATIONS OF INDIVIDUAL SUBSETS 

ZSUM2 = SUM OF SQUARES OF STANDARD DEVIATIONS OF SUBSETS 

NG = NUMBER IN SUBSET GROUP 

NG = NGK 

NR = N/NG 

WRITE (3, 106) NG/NR 

106 FORMAT (//'0 BREAK UP VALUES INTO’ 13 , ' SUBSETS OF '13,' VALUES EACH ) 

NU = NR*NG 

SIZE OF GROUP SUBSET COMPUTATION COMPLETED 

COMPUTE AVERAGE AND STANDARD DEVIATION FOR THE VARIOUS SUBSETS OF VALUES 

SUM = 0 

SUM2 = 0 

GSUM = 0 

GSUM2 = 0 

ZSUM = 0 

ZSUM2 = 0 

NG = 0 

WRITE (3, 110) 

110 FORMAT (' 0 ' 5X, ' SUBSET AVE'21X, 'ITS SIGMA'/' ') 

MAKE STUDY OF SUBSETS 
DO 200 I = 1,K 
SUM = SUM + X(I) 

SUM2 = SUM2 + X(I) **2 
IF (I— NR* (I/NR) ) 200,120,200 
NOTE I (MOD NR) HERE 

COMPLEAT COMPUTATIONS FOR THE GIVEN SUBSET 
120 AVE = SUM/FLOAT (NR) 

SDS = SQRT ( (SUM2— FLOAT (NR) *AVE**2)/ (NR— 1) ) 

GSUM = GSUM + AVE 
ZSUM = ZSUM + SDS 
GSUM2 = GSUM2 + AVE* *2 
ZSUM2 = ZSUM2 + SDS* *2 
WRITE (3, 150) AVE, SDS 
150 FORMATC ' F15 . 4 , 15X, F15 . 4 ) 

NG=NG+1 

HAS FULL SET OF SUBSETS BEEN CONSIDERED 
IF (I— NU) 180,210,180 

ZERO ACCUMULATORS FOR SUMS OF VALUES ASSOCIATED WITH SUBSETS 
180 SUM = 0 
SUM2 = 0 
200 CONTINUE 

210 WRITE (3, 211) mrTT ,-„ „ 

211 FORMAT (///5X* AVE OF AVES THEIR SIGMA AVE OF SIGMAS THExR SI 

.GMA'/' ') 

COMPUTE GRAND AVERAGES AND STANDARD DEVIATIONS 

G MEANS GRAND, Z MEANS STANDARD DEVIATION IN FOLLOWING FOUR NAMES 
GAVE = GSUM/FLOAT (NG) 

ZAVE = ZSUM/FLOAT (NG) 

GAVEZ = SQRT ( (GSUM2-FLOAT (NG) *GAVE**2) / (NG-1) ) 

ZAVEZ = SQRT ( (ZSUM2-FLOAT (NG) *ZAVE**2) / (NG-1) ) 

WRITE (3, 30) GAVE, GAVEZ, ZAVE, ZAVEZ 

RETURN 

END 






c 

c 

c 



c 

c 

c 

c 

c 

c 

c 

c 



PUTS THE TWO COLUMN VECTORS X(I) AND Y(I) IN ORDER ACCORDING TO THE 
VALUES IN X(I) . THE PAIR X(J) AND Y(J) ARE KEPT TOGETHER IN ORDERING. 

THE SMALLEST X(I) IS PUT IN X(l) 

DIMENSION 'X(l) ,Y(1) 

CHECK EVERY ELEMENT 

30 DO 50 J = 1/K 

NOTE LOWER AS INITIAL SMALLEST ELEMENT 

L = J 

CHECK EACH ELEMENT HIGHER THAN THE FIRST 

DO 40 I = J/K 

IF THE LOWER IS ,LE. THE HIGHER ONE ALL O.K. 

IF HIGHER ELEMENT SMALLER THAN LOWER ONE, NOTE THIS 

35 L = I 
40 CONTINUE 



TX = X(J) 

TY = Y (J) 

X(J) = X (L) 
Y(J) = Y (L) 

X(L) = TX 
50 Y (L) = TY 
RETURN 
END 



L IS KEPT THE INDEX OF THE SMALLEST ELEMENT EXAMINED 

ON EXIT FROM DO LOOP L IS INDEX OF SMALLEST ELEMENT 
TEMPORARILY STORE VALUES OF LOWEST ELEMENT STUDIED 

TRANSFER SM ALLE ST ELEMENT TO LOWEST POSITION STUDIED 

PLACE TEMPORARY VALUE IN PLACE OF SMALLEST VALUE 



* 



♦ 
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SUBROUTINE COUNT (L, DX, SP , K, LX, QLB , QUB , AVE ) 

C THIS SUBROUTINE ACCEPTS L VALUES DX(I) . A HISTOGRAM ANALYSIS 

C IS MADE BY SETTING UP K INTERVALS OF WIDTH SP WHICH SPAN THE 

C SET OF VALUES INCLUDED IN DX(I), HAVING UPPER AND LOWER BOUNDS 

C STORED IN QUB (I), AND QLB (I). THE NUMBER OF VALUES DX(I) WHICH ARE 

C FOUND TO FALL IN THE VARIOUS INTERVALS ARE STORED IN LX (I) . 

C A HISTOGRAM OF LX (I) IS INCLUDED USING ************** 

C THE NUMBER OF HISTOGRAM BARS IS AN ODD NUMBER JUST LESS THAN THE 

C SQUARE ROOT OF THE NUMBER OF VALUES IN THE LIST 

INTEGER STAR, BLANK 

DIMENSION DX(1) ,LX(1) ,QLB(1) ,QUB(1) ,ID(30) 

BLANK = 16448 
STAR = 23644 
QL = L 

K = SQRT(QL) - 0.5 
IF (K- (K/2) *2) 4,2,4 
2 K = K - 1 
4 IF (K-3) 6,6,8 
6 K = 3 
8 QK = K 

DO 10 I = 1,K 
10 LX (I) = 0 

SP = 1.0001* (DX(L)-DX(1))/QK 
AM = DX(1) + (DX (L) - DX(l))/2.0 
DO 20 I = 1, K 
QI = I - K/2 - 1 
DO 18 J = 1,L 
QLB (I) = (QI— .5) *SP + AM 
QUB (I) = (QI+.5) *SP + AM 
IF (DX( J) -QLB (I) ) 18,18,14 
14 IF (DX(J) -QUB (I) ) 16,16,18 
16 LX (I) = LX (I) + 1 
18 CONTINUE 
20 CONTINUE 

DO 25 I = 1,30 
25 ID (I) = I 
WRITE (3 31 ) 

31 FORMAT (///‘0*4X, ‘LOWER BOUND ' 4 X, ‘UPPER BOUND HISTOGRAM BARS*/* ‘) 

DO 60 I = 1,K 
LXU = LX (I) 

QLB (I) = QLB (I) + AVE 
QUB (I) = QUB (I) + AVE 
IF (LXU- 90) 37,37,36 

35 FORMAT (‘ ‘ ,2F14.4,90Al) 

36 LXU = 90 

37 WRITE (3, 35) QLB (I) , QUB (I) 

IF (LXU) 55,55,50 

50 WRITE (3, 100) (STAR,L = 1,LXU) 

55 QLB (I) = QLB (I) - AVE 
60 QUB (I) = QUB (I) - AVE 
100 FORMAT (‘+‘,29X,90A1) 

RETURN 

END 
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SUBROUTINE CHI (K, X, EX, DX,DX2 , CHI2 , QLB,QUB ,SX, SEX, SDX, SDX2 , SCHI2 , 
1AVE, S) 

CHI IS USED IN CONJUNCTION WITH SUBROUTINE GAUS. SEE THE LATTER 
FOR THE SIGNIFICANCE OF THE VARIABLES IN CHI, SUBROUTINE CHI 
ACCEPTS K VALUES OF X(I) AND EX (I) , NORMALIZES THE VALUES OF 
EX (I) SO THAT THE SUM OF THEIR VALUES IS EQUAL TO THE SUM OF THE 
VALUES OF X(I), THEN CALCULATES THE VALUES OF SEX, SDX, SDX2 , SCHI2 
DIMENSION X(l) ,DX (1) ,EX(1) ,DX2(1) ,CHI2 (1> ,QLB(1) ,QUB(1) 

SX = 0.0 
SEX =0.0 
DO 10 I = I,K 

SX = SX + X(I) 

10 SEX = SEX + EX (I) 

DO 15 I = 1,K 

15 EX (I) = (SX/SEX) *EX (I) 

SEX =0.0 
SDX =0.0 
SDX2 =0.0 
SCHI2 =0.0 
DO 30 I = 1,K 
SEX = SEX + EX (I) 

DX(I) = EX (I) - X(I) 

IF (ABS (DX(I) )— 1.0E—15) 25,26,26 

25 DX(I) = 1.0E-15 

26 SDX = SDX + ABS (DX(I) ) 

DX2 ( I ) = DX (I) **2 
SDX2 = SDX2 + DX2 (I) 

IF (EX (I) -1 . OE-15) 27,28,28 

27 EX (I) = 1.0E-15 

28 CHI2 (I) = DX2 (I) /EX (I) 

30 SCHI2 = SCHI2 + CHI2(I) 

QK = K 

WRITE (3,100) (X(I) ,EX(I) ,DX(I) ,DX2(K) ,CHI2 (I) ,QLB(I) ,QUB (I) , 1=1, K) 

100 FORMAT (////' NOW LETS LOOK AT THE EXPECTED HISTOGRAM '//' 0 COUNT E 
.XPECT DIF DIF2 CHI2 LOWER UPPER'// (' '7F7.2)) 

WRITE (3, 120) SX, SEX, SDX, SDX2 , SCHI2 
120 FORMAT ( 'O' 7F7.2) 

WRITE (3, 130) 

130 FORMAT ('0 THE LAST LINE OF VALUES GIVES SUMS OF VALUES IN COLUMNS' 
RETURN 
END 
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SUBROUTINE GEX(K,L,EX,SP, S) 

C THERE ARE L VALUES HAVING A GAUSSIAN DISTRIBUTION WHICH HAVE 

C BEEN FOUND TO FALL IN K INTERVALS OF WIDTH SP. THE VALUE OF 

C K IS AN ODD INTEGER. THE STANDARD DEVIATION OF THE L VALUES IS 

C THE NUMBER S. THIS SUBROUTINE CALCULATES THE VALUES OF EX (I) 

C WHICH ARE PROPORTIONAL TO THE NUMBER OF VALUES EXPECTED TO FALL IN 
C EACH INTERVAL OF WIDTH SP. NOTE THAT THE SUM OF THE EX (I) IS NOT 

C EQUAL TO THE VALUE OF L BECAUSE THE EXPECTED GAUSSIAN DISTRIBUTION 

C EXTENDS BEYOND THE RANGE COVERED BY THE K INTERVALS OF WIDTH SP 
DIMENSION EX(1) 

QL = L 

DO 50 I = 1,K 

QI = I - K/2 - 1 

IF (K- (K/2) *2) 15,14,15 

14 QI = QI + 0.5 

15 SUM = 0 

DO 40 J = 1,10 

Q J = J 

QJ = QJ*0 .1 

40 SUM = SUM + EXP (— ( (QI— ,5+QJ— 0 .05) *SP) **2/ (2 ,*S**2) ) 

50 EX (I) = (QL/S/2 . 5066) *SUM* (SP/10 . 0) 

RETURN 

END 



Sample Data Deck for GAUS 
// XEQ GAUS 

DICK ROWE TORSION PENDULUM PART II 

22.2 22.1 22.2 21.9 22.2 21.0 22.0 22.0 22.3 21.9 22.0 22.0 21.8 22.2 

22.4 22.4 22.7 22.2 21.5 22.4 21.8 22.1 21.8 22.0 22.0 22.4 21.8 22.2 
22.1 22.2 21.7 22.2 22.2 22.2 22.3 22.5 21.6 22.6 22.1 22.0 21.9 21.8 

22.0 21.6 22.3 22.5 22.3 22.0 22.3 22.2 9999 



APPENDIX 



Format-Free Input 

When one is working with students of little or no previous 
experience in using a digital computer, one of the major problems 
encountered is that of their properly entering data on punched 
cards. This is true both in the case where the computer program 
is written for the student by the instructor and the case where 
students write their own programs for simple data analysis. Even 
for experienced programmers, entering data on cards according to a 
specified FORMAT can be a serious nuisance. The Coe College FREE 
STYLE input programs are designed to remove this barrier to easy 
use of the digital computer; it can easily make the difference 
between success and total failure in introducing digital computing 
into instructional activities. 

The basic premise of the system is that a student should be 
able to punch into data cards, in proper sequence, the numbers to 
be entered into the computer with practically no further restric- 
tions. In general, any recognizable number is legal. In particu- 
lar: 

1. Numbers may be placed anywhere on a card or cards 

just so long as: a) one or more blanks separate 

different numbers; b) all of each number is on 
just one card; and c) numbers from two different 
"batches" do not appear on the same card. 

2. Numbers are read in "batches" of one or more val- 

ues. Each "batch" corresponds to the execution 
of a statement of the form CALL FREE(X,NVj. This 
causes up to NV values to be read into the loca- 
tions X(l), X (2) , X (3) , . . ., X(NV) . Successive 
cards are searched for numbers until one of the 
following happens : a) NV values have been read 

in — RETURN to calling program; b) 9999 encounter- 
ed on a data card — NV is set equal to the number 
of values previously read in — RETURN to calling 
program — "9999” thus means "end of batch" when NV 
is set big; c) // found in columns one and two — 
monitor trap. 

3. Numbers may contain: a) up to 8 digits; b) alge- 

braic signs may be used freely — the "+" sign is 
optional; c) a decimal point — optional for inte- 
ger values; d) a properly positioned "E" in an 
exponential constant, provided there is no blank 
space before the "E" and it is followed by an ex- 
plicit "+" or "-" sign (for example, 0.317E-8) . 
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4* Illegal characters are treated as blanks, if col- 
umn one of the card is blank . This permits the 
following z a) entries on data cards of the form 
A = 5.6 so the student may remind himself of the 
meaning of the number; b) erasure of keypunch 
errors by overpunching all of the number in gues 
tion with X’s (this greatly speeds up novice use 
of an electric keypunch and permits ready use of 
manual inexpensi e ones in laboratory situations) ; 
c) insertion of pure comment cards and blank cards 
in the data deck for identification and reminder 
purposes. 

5. Comment cards with a "C" in column one are not 
examined for numerical values; instead they are 
printed out on the printer with carriage control 
according to the content of column two • . This 
permits batched data decks with pagination and 
student name header printed at the top of his 
printed output. 

The card reading by these programs is slow but very useful 
with limited amounts of data. Efforts will be made to write fas 
ter versions— both with only FORTRAN and also taking advantage oj. 
the IDEAL subroutines available on the IBM 113 G which permit a 
faster search and overlapped I/O. 

The following pages contain listings of FREE and associated 
programs from the FREE STYLE package, as well as a flow chart for 
FREE, which is the fundamental subroutine. In the flow chart, 
literals such as "+” within a diamond (decision symbol) indicate 
a test as to whether the column under consideration contains that 
symbol. The symbol "b" stands for a blank column. 
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// * MODIF68 12 
// DUP 

4 DELETE FREE 

// FOR 

*ONE WORD INTEGERS 

SUBROUTINE FREE (X,NV) 

DIMENSION X(l) ,NA(80) ,PT(10) 

C 

C FREE STYLE CARD READER 
C 

C CALL FREE (X,NV) WILL READ UP TO NV VALUES INTO X(l) . . . X(NV) 

C THE FIRST ARGUMENT OF THE CALL STATEMENT MUST BE A SINGLY SUBSCRIPTED REAL 

C VARIABLE 

C THE SECOND ARGUMENT OF THE CALL STATEMENT MUST BE AN INTEGER VARIABLE 

C PUNCH VALUES ANYWHERE ON CARD BUT SEPARATE THEM WITH ONE OR MORE BLANKS 

C VALUES MAY HAVE UP TO EIGHT INTEGER DIGITS PUNCHED TO DESIGNATE THEM 

C DECIMAL POINTS ARE HONORED WHERE PUNCHED BUT ARE OPTIONAL 

C VALUES ARE READ OFF CARDS IN ORDER RUNNING FROM LEFT TO RIGHT ON 

C SUCCESSIVE CARDS UNTIL NV VALUES FOUND OR 9999 FOUND 

C YOU MAY PUT AS MANY OR AS FEW VALUES ON EACH CARD AS DESIRED 

C FOUR NINES STANDING ALONE ON CARD SERVE AS SIGN OF END OF VALUES 

C IF 9999 IS ENCOUNTERED NV IS SET TO THE NUMBER OF VALUES PREVIOUSLY READ 

C TO ENTER THE VALUE 9999 INTO THE COMPUTER ENTER IT AS 9999. 

C 

C WITH NV = 4 CALL FREE (X,NV) WOULD READ 1234 OFF NEXT CARD, LOSE 567 

C 1234567 

C WITH NV = 10 CALL FREE (X,NV) WOULD READ 1 2 3 4 5 OFF NEXT CARD AND SET 
C NV = 5 BEFORE RETURN FROM THE CALL 

C 12345 9999 6789 

C 

C LEGAL CHARACTERS ARE 0123456789 .E+- 

C THERE MUST BE NO BLANK BEFORE THE E AND THERE MUST BE A + OR - AFTER 

C E TYPE NUMBERS ARE HONORED IN FOLLOWING FORMS 12.0E+2 4.5E-2 
C LEGAL VALUES 245 +23 -24 -23. -.02 +.03 12E+3 -4.5E-23 1E+.01 45632578 
C ILLEGAL VALUES 2.13.52.45.6 23.4E 3 23.8 E+35 9999 (EXCEPT AS END) 

C I LLE GAL VALUES 23E05 E+03 456328754 45.2547865 +36-25+47+89 0.00000005 

C 

C 

C CARDS HAVING COLUMN ONE BLANK WILL HAVE ALL ILLEGAL CHARACTERS REMOVED 

C BEFORE THE NUMBERS ARE READ OFF IT 

C DATA CARDS WITH ANY PUNCH IN COLUMN ONE EXCEPT C WILL BE READ WITHOUT ANY 

C CHECKING FOR ILLEGAL CHARACTERS 

C THANKS TO ERASURE FEATURE MISPUNCHED CHARACTERS MAY BE 'X-ED' OUT ON CARD 

C THIS MAKES IT PRACTICAL TO USE MANUAL CARD PUNCHES IN THE LABORATORY 

C 

C NOTE THAT COLUMN ONE MUST BE BLANK FOR ERASURE FEATURE TO OPERATE 
C ERASURE OF ILLEGAL CHARACTERS PERMITS DATA CARDS AS THOSE FOLLOWING 

C HERE CONSIDER COLUMN 7 TO BE COLUMN 1 OF A REAL DATA CARD 

C A = 2.3 B = 5.6 C = 4.2 ON FIRST EXPERIMENT 9999 

C FOLLOWING ARE THE TIMES 4 7 9 12 14 16 18 20 9999 

C FOLLOWING ARE THE DISTANCES 5 6 7 9 10 13 45 9999 

C DDATA BY NANCY PHYSICIST ON AUGUST TENTH (DON’T PUT DIGIT FOR DATE) 

C THIS FEATURE WILL HELP NOVICES KEEP TRACK OF WHAT THE VALUES REPRESENT 

C 

C A 'C' IN COLUMN ONE OF A DATA CARD READ BY 'FREE' WILL BE PRINTED ON 
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THE PRINTER BUT OTHERWISE WILL BE IGNORED 

THE FEATURE OF A *C* IN COLUMN ONE vERMxjlS PRINT-OUT TO BE IDENTIFIED 
FOR LABORATORY DATA REDUCTION JOBS HAVING MANY SETS OF DATA FOR ONE RUN 

CARDS ARE READ AT BETWEEN 60 AND 120 CARDS PER MINUTE DEPENDING ON CONTENT 



NCP = 1 

NCP = INCREMENT OF COLUMN COUNTER. SET NCP = 2 FOR ALTERNATE COLUMNS 
NVR = 1 

NVR = NUMBER OF VALUES READ IN 
PT (1) = 1 

PT (I) = POWERS OF TEN = 10** (1-1) 

DO 10 NC = 2,10 
10 PT (NC) = PT (NC-1) *10 
*NOPX = 1 

NOPX = NUMBER OF OPERATION IN MULTIPLICATION OF ENTRY VALUES 
NOPX = 1 MEANS NO OPERATION NOPX = 2 MEANS MULTIPLY BY POWER OF TEN 
20 READ (2, 32) NA 
NANC = NA (1) 

IF (NANC+15552 ) 40,30,40 

DATA CARD WITH *C' IN COLUMN ONE IS PRINTED ON PRINTER 
•C* CARDS HAVE CARRIAGE CONTROL ON COLUMN 2 

USE 'Cl' TO GO TO TOP OF PAGE, 'CO' FOR DOUBLE SPACE, 'C 1 FOR SINGLE 
30 WRITE (3, 32) NA(2),NA 
32 FORMAT ( 8 1A1) 

GO TO 20 

TO CUT STORAGE DEMANDS PUT A GO TO 70 STATEMENT AT 40 AND CUT FOLLOWING 
40 IF(NANC-16448) 70,45,70 
45 DO 68 NC = 2,80 
NANC = NA(NC) 

IF (NANC-16448) 50,68,50 
CHECK TO SEE IF CHARACTER IS A LETTER 

50 IF (NANC+5 824) 56,56,51 
CHECK TO SEE IF CHARACTER IS AN INTEGER 

51 IF (NANC) 68,52,52 
CHECK TO SEE IF CHARACTER IS A *-* 

52 IF(NANC-24640) 53,68,53 
CHECK TO SEE IF THE CHARACTER IS A ' . ' 

53 IF (NANC-19264 ) 54,68,54 
CHECK TO SEE IF THE CHARACTER IS AN IBM 026 '+' 

54 IF (NANC-2 054 4 ) 55,68,55 
CHECK TO SEE IF THE CHARACTER IS AN IBM 029 '+' 

55 IF(NANC-20032) 56,60,56 
CHECK TO SEE IF THE CHARACTER IS AN 'E' 

56 IF (NANC+150 40 ) 61,57,61 

THE 'E' TYPE NUMBERS MUST HAVE A '+' OR FOLLOW THE 'E‘ 

CHECK TO SEE IF THE 'E' IS FOLLOWED BY *-* OR *+* OR EITHER TYPE 
ir 57 IF (NA(NC+1) -20544) 58,68,58 

58 IF(NA(NC+l)-20032) 59,68,59 

59 IF (NA(NC+1) -24640) 61,68,61 

60 NA(NC) = 20544 
GO TO 68 

61 NA(NC) = 16448 

68 CONTINUE 
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70 NC = 0 

C NC = NUMBER OF COLUMN OF CARD BEING CONSIDERED 

80 NC = NC + NCP 

C LOOK FOR FIRST CHARACTER OF NEXT VALUE 

IF (NC-80) 90,90,20 
90 NANC = NA (NC) 

IF (NANC-16448) 100,80,100 
100 ND = 0 

C ND = NUMBER OF DIGITS IN NUMBER CONSIDERED 

NDD = 1 

C NDD = NUMBER OF DECIMAL DIGITS OF NUMBER CONSIDERED 

NDDF = 0 

C NDDF = INCREMENT OF NDD FOR COUNTING NUMBER OF DIGITS BEYOND ' . ' 

NVAL = 0 

C NVAL = INTEGER VALUE OF UP TO FIRST FOUR DIGITS OF NUMBER 

NVAL 2 = 0 

C NVAL2 = INTEGER VALUE OF DIGITS BEYOND FIRST FOUR DIGITS OF NUMBER 

NC2 = 1 

NC2 = NUMBER OF COLUMNS BEYOND FOURTH DIGIT OF NUMBER 
NC2 AND NDD ARE AUGMENTED BY ONE FOR USE IN PT(I) 

NOP XL = NOPX 
NOPX = 1 

NOPXL = VALUE OF NOPX LAST TIME X(I) HAD VALUE STORED IN IT 
NSG = +1 

CHECK FOR ALGEBRAIC SIGN OF VALUE 
IF (NANC-24640) 130,120,130 
120 NSG = -NSG 
GO TO 200 

CHECK FOR BOTH CODES FOR •+• 

130 IF (NANC-20544) 190,200,190 
190 IF (NANC-2 0032) 220,200,220 
200 NC = NC + NCP 

IF (NC-80) 202,202,340 
202 NANC = NA (NC) 

CHECK FOR DECIMAL POINT 
220 IF (NANC-19264) 240,230,240 
230 NDDF = 1 
GO TO 200 
CHECK FOR 'E' 

240 IF (NANC+15040) 260,320,260 

END OF A VALUE ON DATA CARDS INDICATED BY PRESENCE OF A BLANK 
CHECK FOR DELIMITER BLANK 
260 IF (NANC-16448) 270,340,270 

NANC MUST REPRESENT AN INTEGER IF IT REACHES STATEMENT 270 
FIND THE INTEGER VALUE OF THE CHARACTER IN COLUMN NC 
THIS PROCEDURE SPECIFIC TO IBM 1130 CHARAC TER CODE 

ANOTHER ROUTINE HERE MAY BE SUBSTITUTED FOR ANOTHER CHARACTER CODE 
270 I = (NANC+4032) /256 
C COUNT NUMBER OF DIGITS AND DECIMAL DIGITS 

ND = ND + 1 
NDD = NDD + NDDF 
IF (ND-5 ) 275,285,285 
275 NVAL = NVAL* 10 + I 
GO TO 200 

285 NVAL2 = NVAL2*10 + I 
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NC2 = NC2 + I 
GO TO 200 
320 NOPX = 2 

340 X(NVR) = ( (NVAL*PT(NC2) + NVAL2) /PT (NDD) ) *NSG 
NVR = NUMBER OF VALUES READ IN OFF CARD 

IF THERE WAS AN 'E' DELIMITER ON LAST VALUE MULTIPLY BY POWER OF TEN 
FRACTIONAL POWERS PERMITTED HERE 
21.3E+0.5 IS A LEGAL DATA VALUE FOR 'FREE' 

GO TO (360, 350) ,NOPXL 
350 X(NVR-l) = X (NVR-1) *10 . 0**X (NVR) 

NVR = NVR -1 
360 GO TO (410,480) , NOPX 
410 IF(NDDF) 440,430,440 

FOUR NINES IN A ROW WITHOUT DECIMAL POINT INDICATE END OF STRING OF VALUES 
430 IF (NVAL-9999 ) 440,450,440 
READ ONLY UP TO 'NV' VALUES 
440 IF (NVR-NV) 480,460,460 

DO NOT RETURN THE VALUE 9999 
450 X (NVR) = 0 

CALL FREE WITH SECOND ARGUMENT A VARIABLE NEVER A CONSTANT 
NOTE THAT HERE THE SECOND ARGUMENT HAS ITS VALUE CHANGED 
NV ON RETURN IS THE NUMBER OF VALUES ACTUALLY READ IN BEFORE 9999 
NV = NVR - 1 
460 RETURN 
480 NVR = NVR + 1 
GO TO 80 
END 

// DUP 

* STORE WS UA FREE 
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(A Program to Input Parallel Column Vectors) 



SUBROUTINE FREXY (X,Y,NV) 

READS UP TO NV VALUES OF X(I) AND Y(I) IN SEQUENCE END BY 9999 
SEE DOCUMENTATION ON FREE 

IN CALLING PROGRAM X MUST BE DIMENSIONED IN EXCESS OF 2*NV 
IF DIMENSION X(2*K),Y(K) YOU MAY EQUIVALENCE X(K+1) AND Y(l) TO SAVE Ct)RE 
DIMENSION X (2) , Y (1) 

THE DIMENSION STATEMENT MERELY DECLARES THESE TO BE SUBSCRIPTED VARIABLES 
NV = NV*2 
CALL FREE (X,NV) 

NV = NV/2 
IM = NV - 1 

IN CASE ONLY ONE PAIR OF VALUES READ IN 
Y (1) = X (2) 

IF (IM) 500,500,20 

SHIFT EVEN NUMBERED LOCATIONS UP AND ODD NUMBERED LOCATIONS TO BOTTOM 
20 DO 100 I = 1,IM 
12 = 2*1 
I2P =12+1 

FOR EACH EVEN NUMBERED LOCATION SHIFT NEXT ODD ONE OUT TO TEMP FIRST 
THEN SHIFT ALL ORIGINAL EVEN LOCATIONS BELOW IT UP ONE LOCATION 
TEMP = X(I2P) 

DO 80 J = 1,1 
K = I2P - J 
X (K+l) = X (K) 

80 CONTINUE 

C PUT THE ODD LOCATION IN PLACE MADE AT BOTTOM OF EVEN LOCATION STRING 

X(I+1) = TEMP 
100 CONTINUE 

C PUT ALL EVEN LOCATIONS INTO THE Y COLUMN VECTOR 

DO 200 N = 1, NV 
K = NV + N 
200 Y(N) = X (K) 

500 RETURN 
END 
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Subroutines That Complement Free 
By Facilitating Reading In Unscripted Values 



SUBROUTINE RD1 (A) 
DIMENSION X(l) 

NV = 1 

CALL FREE (X,NV) 

A = X(l) 

DpfTITTmT 

IUjXUIUI 

END 



SUBROUTINE IRD1 (NA) 
DIMENSION X(l) 

NV = 1 

CALL FREE (X,NV) 

NA = X(l) +0.0001 

RETURN 

END 



SUBROUTINE RD2 (A,B) 
DIMENSION X (2) 

NV = 2 

CALL FREE (X,NV) 

A = X(l) 

B = X (2) 

RETURN 

END 



SUBROUTINE RD3(A,B,C) 
DIMENSION X (3) 

NV = 3 

CALL FREE (X,NV) 

A = X(l) 

B = X(2) 

C = X (3) 

RETURN 

END 



SUBROUTINE IRD2 (NA,NB) 
DIMENSION X (2) 

NV = 2 

CALL FREE (X,NV) 

NA = X(l) +0.0001 
NB = X (2) +0.0001 
RETURN 
END 



SUBROUTINE IRD3 (NA,NB,NC) 
DIMENSION X (3) 

NV = 3 

CALL FREE (X,NV) 

NA = X (1) +0.0001 
NB = X (2) +0.0001 
NC = X (3) +0.0001 
RETURN 
END 
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(Sample Print-Out From "GAUS" ) 



DICK ROWE TORSION PENDULUM PART II 



OUTPUT FROM SUBROUTINE 


SUBSET 






VALUES 


RUNNING 


RUNNING 


RUNNING 


AS READ IN 


AVERAGE 


SIGMA 


SIGMA OF AVE 


22.2000 


22.2000 






22.1000 


22.1500 


0.0704 


0.0498 


22.2000 


22.1666 


0.0596 


0.0344 


21.9000 


22.0999 


0.1419 


0.0709 


22.2000 


22.1199 


0.1305 


0.0583 


21.8000 


22.0666 


0.1757 


0.0717 


22.0000 


22.0571 


0.1622 


0.0613+ 


22.0000 


22.0499 


0.1513 


0.0535 


22.3000 


22.0777 


0.1646 


0.0548 


21.9000 


22.0599 


0.1648 


0.0521 


22.0000 


22.0545 


0.1575 


0.0475 


22.0000 


22.0499 


0.1514 


0.0437 


21.8000 


22.0307 


0.1604 


0.0445 


22.2000 


22.0428 


0.1611 


0.0430 


22.4000 


22.0666 


0.1809 


0.0467 


22.4000 


22.0874 


0.1938 


0.0484 


22.7000 


22.1235 


0.2395 


0.0580 


22.2000 


22.1277 


0.2334 


0.0550 


21.5000 


22.0947 


0.2686 


0.0616 


22.4000 


22.1099 


0.2703 


0.0604 


21.8000 


22.0952 


0.2720 


0.0593 


22.1000 


22.0954 


0.2653 


0.0565 


21.8000 


22.0825 


0.2665 


0.0555 


22.0000 


22.0791 


0.2611 


0.0533 


22.0000 


22.0759 


0.2561 


0.0512 


22.4000 


22.0884 


0.2592 


0.0508 


21.8000 


22.0777 


0.2605 


0.0501 


22.2000 


22.0821 


0.2569 


0.0485 


22.1000 


22.0827 


0.2520 


0.0468 


22.2000 


22.0866 


0.2488 


0.0454 


21.7000 


22.0741 


0.2543 


0.0456 


22.2000 


22.0781 


0.2511 


0.0444 


22.2000 


22.0817 


0.2485 


0.0432 


22.2000 


22.0852 


0.2454 


0.0420 


22.3000 


22.0914 


0.2449 


0.0413 


22.5000 


22.1027 


0.2505 


0.0417 


21.6000 


22.0891 


0.2608 


0.0428 


22.6000 


22.1026 


0.2701 


0.0438 


22.1000 


22.1025 


0.2662 


0.0426 


22.0000 


22.0999 


0.2633 


0.0416 


21.9000 


22.0950 


0.2620 


0.0409 


21.8000 


22.0880 


0.2630 


0.0405 


22.0000 


22.0860 


0.2603 


0.0396 


21.6000 


22.0749 


0.2672 


0.0402 
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22.3000 
22.5000 

22.3000 

22.0000 

22.3000 

22.2000 

THERE WERE 50 VALUES 



BREAK UP VALUES INTO 

SUBSET AVE 

22.0571 
22.0285 
22.1999 
22.0428 
22.1285 
22.0714 
22.1428 



22.0799 


0.2665 


22.0891 


0.2706 


22.0935 


0.2698 


22.0916 


0.2672 


22.0958 


0.2663 


22.0979 


0.2643 


IN THE FULL 


CLEAN SET 



7 SUBSETS OF 7 VALUES EACH 

ITS SIGMA 

0.1622 
0.1709 
0.4125 
0 ..2152 
0.1979 
0.3640 
0.2992 



0.0397 

0.0399 

0.0393 

0.0385 

0.0380 

0.0373 



0 




AVE OF AVES 
22.0959 



ORDERED VALUES 

21.5000 

21.6000 

21.6000 

21.7000 

21.8000 

21.8000 

21.8000 

21.8000 

21.8000 

21.8000 

21.9000 

21.9000 

21.9000 

22.0000 

22.0000 

22.0000 

22.0000 

22.0000 

22.0000 

22.0000 

22.0000 

22.0000 

22.1000 



THEIR SIGMA 
0.0635 



DEV FROM AVE 

- 0.5979 
-0.4979 
-0.4979 
-0.3979 
- 0.2979 
-0.2979 
-0.2979 
- 0.2979 
-0.2979 
-0.2979 
-0.1979 
-0.1979 
-0.1979 
-0.0979 
-0.0979 
-0.0979 
-0.0979 
-0.0979 
-0.0979 
-0.0979 
-0.0979 
-0.0979 
0.0020 



AVE OF SIGMAS 
0.2603 



THEIR SIGMA 
0.0991 



32 



22.1000 

22.1000 

22.1000 

22.2000 

22.2000 

22.2000 

22.2000 

22.2000 

22.2000 

22.2000 

22.2000 

22.2000 

22.2000 

22.2000 

22.3000 

22.3000 

22.3000 

22.3000 

22.3000 

22.4000 

22.4000 

22.4000 

22.4000 

22.5000 

22.5000 

22.6000 

22.7000 



0.0020 

0.0020 

0.0020 

0.1020 

0.1020 

0.1020 

0.1020 

0.1020 

0.1020 

0.1020 

0.1020 

0.1020 

0.1020 

0.1020 

0.2020 

0.2020 

0.2020 

0.2020 

0.2020 

0.3020 

0.3020 

0.3020 

0.3020 

0.4020 

0.4020 

0.5020 

0.6020 



LOWER BOUND UPPER BOUND HISTOGRAM BARS 



21.4999 

21.7399 

21.9799 

22.2200 

22.4600 



21.7399 **** 

21.9799 ********* 

22.2200 ************************ 

22.4600 ********* 

22.7000 **** 



NOW LETS LOOK AT THE EXPECTED HISTOGRAM 



COUNT 


EXPECT 


DIF 


DIF2 


CHI 2 


LOWER 


UPPER 


4.00 


_ 3.69 


“0.30 


0.09 


0.02 


-0.59 


-0.35 


9.00 


12.21 


3.21 


10.33 


0.84 


-0.35 


-0.11 


24.00 


18.18 


-5.81 


33.75 


1.85 


-0.11 


0.12 


9.00 


12.21 


3.21 


10.33 


0.84 


0.12 


0.36 


4.00 


3.69 


-0.30 


0.09 


0.02 


0.36 


0.60 


50.00 


50.00 


12.85 


54.61 


3.59 






THE LAST 


LINE OF 


' VALUES 


! GIVES 


SUMS OF 


VALUES 


IN COLUMNS 




Flow Chart for FREE 
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INTRODUCTION 



The integration of physics concepts with computer programming 
skills has enabled, students with rudimentary mathematical back- 
ground to solve problems of increasing complexity . Both science 
and non— science students r when required to solve physics problems 
by the application of computer techniques.- appear to grasp, not 
only the basic concepts of physics, but the general capabilities 
and limitations of the computer as well. 

The two sample problems which follow can be solved by students 
who have had only an elementary introduction to analytic geometry 
and differential calculus. The discussion of these problems is in 
two parts: part one deals with the description of sinple harmonic 

motion; part two deals with damped simple harmonic motion. Experi- 
mentally, students in the second year of the science curriculum at 
the Naval Academy were found to spend two to three terminal hours 
on this material after having had one hour of experience with the 
BASIC programming language and an introductory lecture on the phys- 
ics involved. 



STUDENT MANUAL 



Simple Harmonic Motion 

Any motion that repeats itself in equal intervals of time is 
called periodic, or harmonic. If a body moves about an equilibri- 
um position due to a force that is proportional to the distance 
from the equilibrium position to the body itself, then the body is 
said to undergo simple harmonic motion. An ideal example of such 
a system is a block set on a frictionless plane and attached to a 
spring. (See Figure 1.) The force on the block due to the spring 
is always such as to pull the block back to its equilibrium posi- 
tion, (X = 0) and is properly called the restoring force. At 
equilibrium, of course, the force of the spring on the block is 
equal to zero. At any particular instant, the restoring force is 

F = -K*X 

where the asterisk (*) indicates multiplication. The minus sign 
in the equation shows that this force will always be directed op- 
posite to the displacement in the X-direction. 

If Newton's Second Law, F = M*A, is applied to the motion of 
the block of mass M, and if the force F is replaced by the ex- 
pression -K*X, the following relationship results: 

— K*X = M*A = M*(d 2 X/dT 2 ) 

since A (acceleration) is the second derivative of X with re- 
spect to time, T. Thus 

M*(d 2 X/dT 2 ) + K*X = 0 (1.1) 

which describes the periodic motion of the block-spring system. 

From this differential equation of motion, the problem is to 
determine the position of the mass at every instant of time after 
the system has been given an initial displacement. Equation 1.1 
is a second-order differential equation for which a numerical so- 
lution can be obtained if a set of initial conditions are known. 
For instance, in this example, it will be assumed that at time 
T = 0, the instantaneous velocity of the block (dX/dT) is zero, 
and the displacement x = 25 meters. (For idealized experiments, 
the cost of equipment is no object.) 

The first step in the solution of this problem is to re-write 
Equation 1.1 as a system of first-order differential equations. 

The solution of the first-order equations is relatively simple as 
long as one remembers that the expression dX/dT is the s lope of 
the X versus T curve at any instant of time. To re-write the 
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above second-order equations, recall the following definitions of 
instantaneous velocity and acceleration: 

V = dX/dT ; A = dV/dT = d 2 X/dT 2 



It is then possible to re-write Equation 1.1 as follows: 

M* (dV/dT) + K*X = 0 

or dV/dT = -K*X/M 

and dX/dT = V 



( 1 . 2 ) 



(1.3) 



The problem can now be solved approximately if the incremental 
changes in velocity (AV) and displacement (AX) which occur during 
small changes in time (AT) can be approximated. In other words, 
after a small change in time, (AT in BASIC will be called D) , the 
new displacement and velocity can be found by 



X = X + AX 

V = V + AV 

A = — K*X/M 

These equations are in the form of "assignment statements, X — 

X + AX, which direct the computer to calculate a new value o 
by adding AX to the current value. 

A crude first-order method for approximating the change in 
displacement and velocity would be 

AV = A*D 

(change in velocity) = (acceleration) * (change in time) 

AX = V*D 

(change in displacement) = (velocity) ‘(change in time) 

where A and V represent the acceleration and velocity respec- 
tively at the beginning of the time interval. Therefore, 

X = X + V*D 

V = V + A*D 
A = -K*X/M 



A much better estimate could be made by using a weighted av- 
eraae of A and V which could be obtained from values at the 
beginning, midpoint, and end of the time interval. (See Figures 
2a and 2b.) In this figure, the initial slopes 



XI = X 
VI ~ V 
Al — — K*Xl/M 
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are used to predict the slopes, V2 and A2, at the midpoint of the 
time interval, AT = D. Curved arrows indicate slopes of tangents 
to the X— T and V— T graphs, 

X2 = XI + Vl*D/2 
V2 = VI + Al*D/2 
A2 = -K*X2/M 

Next, the values of V2, A2 are used to make a second predic- 
tion of the midpoint values of the slopes V3 and A3 itarting 
from XI and VI. (See Figures 3a and 3b.) 

X3 = XI + V2*D/2 
V3 = VI + A2*D/2 
A3 = -K*X3/M 

Finally, this set of slopes is used to predict the values of the 
slopes V4 and A4 at the end of the time interval (D) , starting 
from XI and Vl. (See Figures 4a and 4b.) 

X4 = XI + V3*D 
V4 = VI + A3*D 
A4 = -K*X4/M 

The theory of the Runge-Kutta approximation shows that if the 
slopes are weighted by factors 1/6, 2/6, 2/6 and 1/6, respectively, 
then the approximation will be fourth— order; i.e., correct to within 
errors proportional to D 5 . Thus, if we set 

V = (VI + 2*V2 + 2*V3 + V4)/6 

A = (Al + 2*A2 + 2*A3 + A4)/6 

then the new values of displacement and velocity are found by using 
the weighted averages as follows: 

X = X + V*D = X + (Vl + 2*V2 -f 2*V3 + V4)*D/6 

V = V + A*D = V + (Al + 2*A2 + 2*A3 + A4)*D/6 

Note that the values of XI, X2, X3 and X4 were computed in order 
to obtain the corresponding slopes Al, A2, A3 and A4 of the V-T 
curve, which in turn predicts the next value of V which predicts 
the next value of X and so on. 

The above algorithms have been incorporated into the computer 
program SIMPLE HARMONIC MOTION for the determination of the displace- 
ment, velocity and acceleration at equal increments of time after the 
mass is set in motion. This program is designed to give graphic out- 
put of X(T) directly, however, the PRINT statement, number 300, 
can be changed to 



300 PRINT T, X, V 
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T*D/2 | I T+0/2 i j T+D/2 | 

T*D T T+D T TfD 
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for m^nerical output. TO, XO, and VO are initial values of T, 
X and V. 



100 REM *************** SIMPLE HARMONIC MOTION *************** 
120 REM 

130 REM *************** INPUT DATA 
140 REM 

150 READ K, D, TO, X0, V0, M 

160 DATA 2.8, 0.1, 0, 25, 0, 1 

170 PRINT "SPRING CONSTANT = " K"N£WTONS /METER" 

180 REM 

190 REM ************ INITIAL CONDITIONS 

200 REM 

210 LET T = TO 

220 LET X = X0 

230 LET V = V0 

240 LET N = 0 

250 REM 

260 REM *************** PLOT ROUTINE 
270 REM 

280 PRINT " X - DISPLACEMENT IN METERS" 

290 PRINT- 

300 PRINT TAB (33); "0" 

310 PRINT "SECONDS 

320 PRINT " " 

330 LET Y1 = INT (X + .5) +36 
340 PRINT T; TAB (Yl) ; "*" 

350 LET N = N + 1 
360 IF N = 40 THEN 580 
370 REM 

380 REM ************* CALCULATE SLOPES 

390 LET XI = X 

400 LET VI = V 

410 LET A1 = — K*X1/M 

420 LET X2 = XI + Vl*D/2 

430 LET V2 = VI + Al*D/2 

440 LET A2 = -K*X2/M 

450 LET X3 = XI + V2*D/2 

460 LET V3 = VI + A2*D/2 

470 LET A3 = -K*X3/M 

480 LET X4 = XI + V3*D 

490 LET V4 = XI + A3*D 

500 LET A4 = -K*X4/M 

510 REM ********** NEW VALUES OF X, V, AND A 

520 LET X = X + (VI + 2*V2 + 2*V3 + V4)*D/6 

530 LET V = V + (A1 + 2*A2 + 2*A3 + A4)*D/6 

540 LET A = -K*X/M 

550 LET T = T + D 

560 SETDIGITS (3) 

570 GO TO 330 
580 END 
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SPRING 


CONSTANT = 2.8 


NEWTONS/METER 


TIME 


X-DISPLACEMENT X-VELOCITY 


(SECS.] 


1 (METERS) 


(METERS/SECOND) 


0 


30 


0 


.05 


29.9 


-4.2 


.1 


29.6 


-8.36 


.15 


29.1 


-12.5 


.2 


28.3 


-16.5 


.25 


27.4 


-20.4 


.3 


26.3 


-24.2 


.35 


25. 


- 27.7 


.4 


23.5 


-31.1 


.45 


21.9 


-34.3 


.5 


20.1 


-37.3 


.55 


18.2 


-39.9 


.6 


16.1 


-42.3 


.65 


13.9 


-44.5 


.7 


11.7 


-46.2 


.75 


9.32 


-47.7 


.8 


6.9 


-48.9 


.85 


4.44 


-49.6 


.9 


1.94 


-50.1 


.95 


-.566 


-50.2 


1. 


-3.07 


-49.9 


1.05 


-5.55 


-49.3 


1.1 


-8. 


-48.4 


1.15 


-10.4 


-47.1 


1.2 


-12.7 


-45.5 


1.25 


-14.9 


-43.5 


1.3 


-17.1 


-41.3 


1.35 


-19.1 


-38.8 


1.4 


-20.9 


-36. 


1.45 


-22.6 


-32.9 
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SPRING CONSTANT = 2.8 NEWTONS/METER 

X - DISPLACEMENT IN METERS 



SECONDS 

0 



0 



* 

* 



1 . 

1.1 

1.2 

1.3 

1.4 

1.5 

1.6 

1.7 

1.8 

1.9 
2 . 
2.1 
2.2 

2.3 

2.4 

2.5 

2.6 

2.7 

2.8 

2.9 
3. 
3c 1 

3.2 

3.3 

3.4 

3.5 

3.6 

3.7 

3.8 

3.9 
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Damped Harmonic Motion 

In a real spring problem, such as one may study in the labora- 
tory, the amplitude of oscillation gradually decreases to zero. 
This is, of course, due to friction, and can be accounted for by 
the addition of damping forces. The oscillations that result are 
called damped harmonic motion. In Figure 5 friction has been sim- 
ulated by taking a disk, which is attached to the spring-mass sys- 
tem, and immersing it in a fluid. As a first approximation, one 
could assume that the faster the disk moves, the greater will be 
the force of friction, or, in equation form: 

friction force = - (constant) * (velocity in X-direction) 



or, F = -B* (dX/dT) 

The minus sign indicates that the friction force will be in a di- 
rection opposite to the direction of motion at any instant of time. 

In the free-body diagram, it is clear that both the spring 
force (-K*X) and the friction force (-B* (dX/dT)) act upon the mass, 
M, at all times to produce the resulting damped harmonic motion. 
Again, beginning with Newton’s Second Law, 

F = M*A 

— K*X - B* (dX/dT) = M*(d 2 X/dT 2 ) 

dividing through this equation by M, and rearranging terms, yields 

d 2 X/dT 2 + B/M* (dX/dT) + K/M*X =0 (2.1) 



(Note: the computer reads K/M*X as (K/M)*X.) Equation 2.1 is a 

second-order differential equation which can be very effectively 
solved by numerical methods if proper initial conditions are given. 
As in the previous problem. Equation 2.1 can be transformed to a 
combination of first-order differential equations. First, by def- 
inition, dX/dT = V, velocity, and dV/dT = d 2 X/dT 2 = A, so that 

dV/dT = — B/M*V - K/M*X (2.2) 

In this manner, the problem has been reduced to one which contains 
two first-order differential equations which can be solved, as be- 
fore, by a series of approximations based on initial conditions and 
estimates of the slopes at each increment of time AT (in BASIC AT 
will be replaced by D) . 

Therefore, after each increment of time. 



and 



X = X + (VI + 2*V2 + 2*V3 + V4)*D/6 
V = V + (A1 + 2*A2 + 2*A3 + A4)*D/6 



and, from Equation 2.2, 



(2.3) 



(A1 + 2*A2 + 2*A3 + 
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Figure 5 
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A = -B/M*V - K/M*X 

The slopes, VI through V4 and Al through A4 (velocities and 
accelerations, respectively) at the four points shown in Figures 
2, 3, and 4, are found as follows: 

(let Cl = -B/M and C2 = -K/M) 

XI * X 

Vl = V 

Al = C1*V1 + C2*X1 (from Equation 2.2) 

These slopes predict the slopes at the midpoint of the time incre- 
ment, D, by 

X2 = XI + Vl*D/2 

V2 = VI + Al*D/2 

A2 = C1*V2 + C2*X2 

A second estimate of the slopes at the midpoint is made by 

X3 = XI + V2*D/2 

V3 * VI + A2*D/2 

A3 * C1*V3 + C2*X3 

and, finally, the slopes at the end of the time interval are found 
by 

X4 * XI + V3*D 

V4 = VI + A3*B 

A4 = C1*V4 + C2*X4 

The only major difference between this and the undamped case is in 
the calculation of A. Substitution into Equation 2.3 then yields 
the desired values of X and V. 

The DAMPED HABMONIC MOTION program which follows illustrates 
the use of these algorithms to describe the motion of the system. 
For the sake of clarity, the asterisks in the graphic output have 
been connected by hand-drawn straight lines. You can easily com- 
pare the results for different choices of parameters by superimpos 
ing one graph over another on a light-box or against a windowpane. 
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REM *************** DAMPED harmonic motion *************** 
REM 

BEM *************** INPUT DATA 
REM 

READ K, B, M, D, VO, TO, XO 
DATA 2.8, 0.4, 5, 0.10, 0, 0, 30 
PRINT 

PRINT "SPRING CONSTANT = "K"NEWTONS/METER" 

PRINT 

PRINT "DAMPING COEFFICIENT = "B" NEWTONS/METER/SEC" 

REM 

EEjj ************ INITIAL CONDITIONS 
REM 



LET Cl = -B/M 
LET C2 = -K/M 
LET X = X0 
LET V = V0 

let t = to 



*********** PRINT COMPUTER OUTPUT 



REM 
PRINT 

PRINT "TIME" , " X- DI SPLACEMENT " , 
LET M = 0 
PRINT T, X, , V 
IF (T - 30) > = 0 
REM ***** 

LET M = M + 1 



THEN 570 
* * * 



'■ VELOCITY" 



***** CALCULATE SLOPES 



LET XI 
LET VI 
LET Al 
LET X2 
LET V2 
LET A2 
LET X3 
LET V3 
LET A3 
LET X4 
LET V4 
LET A4 
REM 
REM 
REM 
LET 
LET 
LET 
LET 



X 

V 

C1*V1 + C2*Xl 

XI + Vl*D/2 

VI + Al*D/2 
Cl*V2 + C2*X2 
XI + V2*D/2 
VI + A2*D/2 
C1*V3 + C2*X3 
XI + V3*D 
Al + A3*D 
C1*V4 + C2*X4 

********* NEW VALUES OF X, V, A, AND T 



X 

V 

A 

T 



X + (VI + 

V + (Al + 
C1*V + C2*X 
T + D 



2*V2 + 2*V3 + V4)*D/6 
2*A2 + 2*A3 + A4) *D/6 



[F M = 10 THEN 310 



GO TO 330 
END 
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SPRING CONSTANT =2.8 NEWTONS/METER 
DAMPING COEFFICIENT = .4 NEWTONS/METER/SEC 



TIME 


X-DISPLACEMENT 


VELOCITY 


0 


30 


0 


1 . 


24.9499 


-10.1607 


2. 


11.5866 


-15.9949 


3. 


-4.7886 


-15.801 


4. 


-18.2324 


-10.1108 


5. 


-24.2815 


-1.33245 


6. 


-21.3957 


7.23454 


7. 


-11.2697 


12.6184 


8. 


2.00713 


13.1865 


9. 


13.5613 


9.11154 


10. 


19.4955 


2.17251 


11. 


18.173 


-4.9898 


12. 


10.6138 


-9.86009 


13. 


-6. 50444E-2 


-10.9162 


14. 


-9.89872 


-8.08357 


15. 


-15.5225 


-2.64969 


16. 


-15.2991 


3.28984 


17. 


-9.75677 


7.62444 


18. 


-1.23836 


8.9659 


19. 


7.05587 


7.07686 


20. 


12.2503 


2.86502 


21. 


12.7719 


-2.02169 


22. 


8.79869 


-5.82688 


23. 


2 . 06266 


- 7.30666 


24. 


-4.87397 


-6.12401 


25. 


- 9,57636 


-2.89649 


26. 


-10.5765 


1.09269 


27. 


-7.81064 


4.39351 


28. 


-2.5336 


5.9077 


29. 


3.22067 


5.24474 


30. 


7.40841 


2.80357 



O 
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SPRING CONSTANT =2.8 NEWTONS/METER 
DAMPING COEFFICIENT = .4 NEWTONS/METER/SEC 

X - DISPLACEMENT IN METERS 

0 





* 
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SPRING CONSTANT = 5.6 NEWTONS/METER 
DAMPING COEFFICIENT = .4 NEWTONS/METER/SEC 

X - DISPLACEMENT IN METERS 



0 
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SPRING CONSTANT =5.6 NEWTONS/METER 
DAMPING COEFFICIENT = .8 NEWTONS/METER/ SEC 

X - DISPLACEMENT IN METERS 
0 



SECONDS 






TEACHER'S GUIDE 



In this presentation, the functional relationships have been 
developed in BASIC in order to facilitate the relation of the al- 
gorithms describing the physics directly to the computer program. 
For instance, in the first program initial conditions are defined 
in line numbers 210 through 240 of the program. The plotting of 
computer-generated output is directed by line numbers 280 and 340. 
The four slopes estimating Velocity and acceleration are calcula- 
ted in line numbers 390 through 500. Time is incremented in line 
550. The statement SETDIGITS (3) in line 560 is used to set the 
size of the numbers to be printed by the computer, in this case, 
three significant digits. It is included solely for esthetic rea- 
sons. 

Students should be encouraged to write a similar program and 
to explore the effects on the period of changing the spring con- 
stant, K, the mass, M, the initial displacement, X0, or the initial 
velocity, V0. New values of K, M, X0, and V0 can be entered in 
the data statement in line 160. The period can actually be deter- 
mined from the plotted computer output, thus, allowing students to 
graph the period versus any of the above varia bles to determine an 
empirical relationship. (Recall that T = 2ir/M/K . ) 

In the second program numerical output is provided. However, 
to enter the plot routine it is only necessary that line numbers 
280 through 340 of the previous program must be added. Compare the 
numerical output of this program (numbers are printed to six digits) 
to the previous program, where the statement SETDIGITS (3) was used. 
Students should be encouraged to investigate the effect on period 
and amplitude of manipulating values of the. spring constant, mass 
and damping coefficient. 

Solution of the differential equation describing the motion 
of the damped spring-mass system can be found in most handbooks. 

For comparison of results, students should program both the numer- 
ical solution and the closed- form solution of this problem, which, 
for TO = 0, is 

X = e " BT ^ M [X0 coswT + C sinuT] 
where w = /k7m, C = VO/w + (B/M) (XO/u) 

For the sake of realism, these equations have not been scaled 
to dimensionless variables. However, programming could have been 
simplified by scaling time to units of one period, T = t/T. 
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TWO EXPERIMENTS 



Conservation of Momentum 
and 

Simple Harmonic Motion 



David T. Grimsrud 

Muhlenberg College 
Allentown, Pennsylvania 




INTRODUCTION 



The two programs discussed here illustrate the use of the com- 
puter as a computational adjunct to an introductory physics demon- 
stration or laboratory. The programs are intended" to aid the stu- 
dent in analyzing the data acquired in the course of the experiment. 
y ie f irs t program, MOMEN, was developed to give immediate answers 
for a lecture demonstration of the conservation of momentum. This 
demonstration will be part of a one-semester "Physics for Poets" 

which does not have an associated laboratory. The students, 
for the most part, are second-year students who have no preparation 

m physics, but may~ have studied another laboratorv science for one 
year. 



-Since this program has not been tested with any groups of stu— 
dents, we omit mention of any student materials, except to say that 
the students, by this nuncture, will have been introduced to the 
computer and will understand simple programming concepts. The 

of conservation of momentum will have been discussed pri- 
or to the demonstration. P 



The second program, PEND1, was written for use in the general 
physics laboratory, in which single harmonic motion is studied by 
recording the position of a pendulum bob (see Figure 1) as it moves 
through one cycle of its oscillation. This program has been used 
Y approximately 60 first— and third— year science majors concurrent- 
ly studying calculus. The method of calculating velocities and 
accelerations used in the program had been used previously in the 
course, although the students had no previous experience with the 
computer. The program was introduced in a written description giv- 
en to the class a week prior to the experiment. The three-hour 
laboratory period was devoted to a summary explanation of the ex- 
periment and program, obtaining the data, keypunching these data, 
running the program, and plotting the results. 

Inasmuch as these two programs have had little classroom test- 
ing, their student syllabi are not fully developed. However, a 
student manual for the simple harmonic- oscillator is presented in 
the Appendix to this paper. 
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TEACHER* S GUIDE 



Conservation of Momentum 

Experiments in the conservation of momentum in one dimension 
are quite common in the classroom and can take many forms: a linear 
air track may support two colliding gliders, two steel balls may be 
hung on bifilar suspensions, etc. In such an experiment some device 
must be provided to measure the position of the particles at equal 
intervals of time. An open- shuttered Polaroid camera and strobe il- 
lumination, or the spark- trace attachment found on most air tracks, 
would be appropriate. The positions indicated by these devices 
permit calculation of velocities and thence momenta? the momentum 
before collision is compared with the momentum after collision. 

The program is divided into three parts: the velocity of each 
particle before and after collision is calculated? the momenta are 
determined? and the momentum of the system before the collision is 
compared with the momentum after the collision. 

The input data required by the program are: 

1. Masses of two particles, M(l) and M(2) . 

2. The time interval, DELT, between successive position 
measurements . 

3. The number, L, of position measurements before colli- 
sion (chosen such that the number after collision is 
also L) . 

4. The positions of the particles which form the three- 
dimensional array DISP(I, J, K) . 

a. DISP (I, 1, 1) 1=1, . . . , L give the positions 

of particle 1 before collision. 

b. DISP (I, 2,1) 1=1, . . . , L give the positions 

of particle 2 before collision. 

c. DISP (I, 1,2) 1=1, . . . , L give the positions 

of particle 1 after collision. 

d. DISP (1, 2, 2) 1 = 1, . . . , L give the positions 

of particle 2 after collision. 

The technique used to compute the velocity eliminates a common 
error the student might make if he were performing the calculation 
without direction. He is faced with the problem of finding the av- 
erage interval between displacements separated by equal intervals 
of time. If we label these positions X(l), X( 2), . . ., X(K), 

. . . and the intervals D(l) = X(2)-X(l)? D(2) = X(3)-X(2)? D(I) = 
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X (1+1)— X (I) , his intuition will likely direct him to consider 

D = i | D(I) 

11 1=1 

as the average value of the interval. Substitution of the measured 
points X (1) , . . ., X(n+1) into the above expression shows that it 
reduces to 



D = i ( X(n+1) -X (1) ) 

i which neglects most of the data obtained. 

i 

* 

j A more efficient technique (in the sense that each data point 

! is used once to find the average interval) is effected if the posi- 

| tion data are divided into two groups and paired to produce inters 

| vals roughly one-half the total interval. These half-intervals are 

j then averaged, the time required to travel a half-interval is com- 

j puted, and the average velocity is found. 

j Assume, as an example, that nine position measurements (L = 9) 

are made before the collision X(l), X(2), . . ., X(9). The pair- 
i ings chosen to evaluate the average half-interval in this case would 

be X (9) -X(4) , X (8) -X (3) , X(7)-X(2), and X(6)-X(l). To choose the 
proper position index to subtract from X(9) we use integer division, 
dividing L by 2 to give the integer variable L2. In this ex- 
ample, L divided by 2 gives 9/2 which is equal to 4 in integer 
i arithmetic. The sum of the four half-intervals for the first parti- 

cle ( X(9)-X(4) ) + ( X (8) -X (3) ) + ( X (7) -X (2) ) + ( X(6)-X(l) ) 
is computed and stored as the variable named INT0T(1,1). The aver- 
age of these four determinations is called AVINT(1,1) and is found 
by dividing by the real value of L2. Finally, since the value of 
A VINT (1,1) contains the distance traveled in five time intervals, 
we find the velocity of the first particle by dividing AVINT(1,1) 

| by the number of time intervals multiplied by the time per interval. 

! 

It would have been preferable to perform data reduction in 
i terms of standard deviations, thereby utilizing fully the data col- 

lected,* however, the sophistication of the students in this course 
does not admit of it, hence the above stratagem, which is clearly 
preferable to being apprehended in the act of taking useless data. 

The remainder of the computation proceeds very simply after veloci- 
ties have been determined; momenta are calculated and values of 
momenta before and after collision are compared. 

MOMEN and a sample of its output are shown on the following 
pages . 



*"Data Reduction," Paul A. Smith, in Computer - Based Physics : 

An Anthology , R. Blum, et al, published by the Commission on College 
Physics (1969) . 
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Progr am : MOMEN 



C 

C 

C 

C 

C 



90 

80 

100 



= INTOT (J,X) + DISP (L, J,K) - DISP(L2,J,K) 
100,100,80 



T’*E 3-DIMENSIONAL ARRAY DISP(I,J,K) CONTAINS INFORMATION ABOUT THE 
DISPLACEMENTS OF THE TWO MASSES INDEXED BY J . K IS ONE FOR DIS- 
PLACEMENTS BEFORE COLLISION, TWO FOR THOSE AFTER. MASSES ARE DE- 
NOTED BY M(J). L GIVES THE NUMBER OF APPROXIMATELY EQUALLY SPACED 
DISPLACEMENTS SEPARATED IN TIME BY THE INTERVAL DELT 
REAL M(2) , INTOT (2,2), MOMl, MOM2 
DIMENSION DISP (15,2,2) , AVINT(2,2), VEL(2,2) 

READ (2,10) M, DELT, L 
10 FORMAT (2F5 . 1, F4.3, 12) 

READ (2, 12) ( ( (DISP (I, J,K) , 1= 1,L), J= 1,2),K= 1,2) 

12 FORMAT (2F4.1) 

USE INTEGER DIVISION TO DP/IDE THE INPUT DATA INTO TWO PARTS 

L2 = L/2 
RL2 = L2 

FIX VALUES OF L AND L2 SINCE THEY WILL CHANGE IN THE DO LOOPS 
N = L 

— jj2 

COMPUTE TOTAL HALF INTERVAL VALUES 
DO 100 J = 1,2 
DO 100 K = 1,2 
INTOT (J,K) =0.0 
L = N 
L2 = N2 
INTOT (J,K) 

IF (L2 - 1) 

L = L-l 
L2 = L2 - 1 
GO TO 90 
CONTINUE 

COMPUTE AVERAGE VELOCITIES 
TIME = DELT* ( (N+l)/2) 

DO 200 J = 1,2 
DO 200 K = 1,2 
AVINT (J,K) = INTOT (J,K)/RL2 
VEL (J,K) = AVINT (J,K) /TIME 

FORMAT ( ' VELOCITY OF THE FIRST PARTICLE BEFORE COLLISION WAS , 
XF6.1,/5X, ' WHILE THAT OF THE SECOND PARTICLE WAS’,F6.1) 

format^ ' Velocities' of the particles after collision were', 

XF6.1,' AND' ,F6 • 1* RESPECTIVELY.') 

WRITE (3,16) VEL(1,2) , VEL(2,2) 

18 fo^at " ( totIl L moAentum M of > thI^ystem before COLLISION WAS', 
XF8.1' 

WRITE (3,18) MOMl 

20 SSEL; after collision, on the 

XOTHER HAND, WAS ', F8.1) 

WRITE (3,20) MOM2 

DISCR = ((MOMl - MOM2) / (MOMl + MOM2) ) * 100. 



200 

14 



16 
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22 FORMAT (' 1’HUS OUR SIMPLE EXPERIMENT HAS SHOWN THAT MOMENTUM IS CON 
XSERVED WITHIN ' , F6.1) 

WRITE (3,22) DISCR 
WRITE (3, 24) 

24 FORMAT (' PERCENT. WHAT EFFECTS CONTRIBUTE TO THIS DISCREPANCY.') 
CALL EXIT 
END 



Output: 



VELOCITY OF THE FIRST PARTICLE BEFORE COLLISION WAS 9.9 
WHILE THAT OF THE SECOND PARTICLE WAS 0.0 
VELOCITIES OF THE PARTICLES AFTER COLLISION WERE 3.3 AND 13.3 RESPECTIVELY. 
TOTAL MOMENTUM OF THE SYSTEM BEFORE COLLISION WAS 997.7 

TOTAL MOMENTUM OF THE SYSTEM AFTER COLLISION, ON THE OTHER HAND, WAS 1007.7 
THUS OUR SIMPLE EXPERIMENT HAS SHOWN THAT MOMENTUM IS CONSERVED WITHIN -0.4 
PERCENT. WHAT EFFECTS CONTRIBUTE TO THIS DISCREPANCY. 



Simple Harmonic Motion 

This program was written for use in the general physics labor- 
atory at Muhlenberg. Simple harmonic motion is studied by making 
a record of the position of a pendulum bob as the bob moves through 
one cycle of its oscillation. A sketch of the apparatus is shown 
in Figure 1. A spark- trace record of the position of the bob is 
left on sensitized paper placed beneath the pendulum on a curved 
support. The student measures positions from the paper, plots 
these results, computes velocities from displacements, plots these 
values, and finally computes and plots accelerations. The three 
graphs obtained are compared to the trigonometric functions which 
describe the system. 




Figure 1 



The input for the program is: 

1. The time interval in seconds between successive sparks, 

DELT. 

2. The total number of data points, J. 

3. The ordered positions, in centimeters, of the J points. 

We require that the position of the J points be measured from 
the midpoint of the motion (which was marked on the sensitized paper) 
so that the data have the approximate form r sin («t) where R = 
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amplitude of oscillation. A typical trace is sketched in Figure 2. 
Furthermore, the initial data point should be the last negative val- 
ue of the position (1) of the bob before it swings positive. The 
J position coordinates are punched on cards using a free style in- 
put, i.e., the actual number (including its decimal point in the 
appropriate place) is punched, a space, another number, space, etc.* 

The initial computation in the program is the angular frequency 
of the motion. The first zero of displacement from the vertical 
occurs (by design) between the first two points (1) and (2) . This 
location is found by linearly interpolating between these points. 
Since the points then become positive, the next zero occurs just be- 
fore the next subsequent negative position (3) . Assuming that the 
midpoint of the motion has been correctly identified, the time be- 
tween the two zeroes is one-half the period. 




Figure 2 

The phase angle and approximate amplitude R are then calcu- 
lated, the latter simply being identified as the maximum of the 
absolute values of the displacement. Differences between successive 
positions are determined and, by multiplying by the reciprocal of 
the time interval, the average velocity in that time interval is 
found. Similarly, differences in velocity multiplied by the recip- 
rocal of the time interval give average accelerations. 

For purposes of plotting graphs the expected values of the 
displacement R sin (wt) , velocity Rw cos (wt) , and acceleration 
-Rw* sin(wt) are computed at times corresponding to the times for 
which the experimental values have been calculated. 

The six columns of print-out are double spaced. In this way 
the velocities which approximate the instantaneous velocity at the 
midpoint of the interval between two displacements cam be placed on 
lines which alternate with the displacements. In a similar fashion 
the accelerations are printed on lines which alternate with the ve- 
locities . 



ERIC 



*Smith, Paul A., ojd. cit . 
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Program : PEND 1 



C 



C 

C 



C 



C 



C 



C 

C 



DIMENSION X(90) , V(90), DIFF(90), XSIN(90), DIF2(90), A(90),ASIN(9 
10), VCOS (90) 

READ (2,10) DELT, J 
10 FORMAT (F5.3, 12) 

IF (DELT) 50,50,12 
12 IF (J) 50,50,14 
14 CALL FREES (X,J) 

CALCULATE THE ANGULAR FREQUENCY OF THE OSCILLATION 
ZERO = X(1)/(X(1) - X(2) ) 

J2M = (J/2) - 4 
DO 60 I = J2M, J 
IF (X(I) ) 70,70,60 
60 CONTINUE 

70 ZEROl = I - 2 + X(I-1)/(X(I-1) - X(I)) 

T2 = (ZEROl - ZERO ) * DELT 
OMEGA = 3. 141593/T2 

CALCULATE THE PHASE ANGLE OF THE MOTION 
PHI = ZERO * OMEGA * DELT 

CALCULATE THE APPROXIMATE AMPLITUDE OF THE MOTION 
XMAX = X (1) 

DO 100 I = 2, J 

IF (ABS (X (I) ) - XMAX) 100,100,90 
90 XMAX = ABS (X(I) ) 

100 CONTINUE 

CALCULATE THE AVERAGE VELOCITY BETWEEN X(I) AND X(I+1) 

J1 = J - 1 
TDELT = 1.0 / DELT 
DO 200 I = 1, J1 
DIFF(I) = X(I+1) - X(I) 

V(I) = DIFF(I) * TDELT 
200 CONTINUE 

CALCULATE T HE AVERAGE ACCELERATION BETWEEN V(I) AND V(I+1) 

J2 = J - 2 

DO 300 I = 1, J2 

DIF2 (I) = V(I+1) - V(I) 

A (1+1) = DIF2 (I) * TDELT 
300 CONTINUE 

COMPUTE THE EXPECTED VALUES OF DISPLACEMENT AND ACCELERATION 

ARG = OMEGA * DELT 

OMEG2 = OMEGA * OMEGA 

DO' 400 I = 1,J 

RLI = I 

XSIN(I) = XMAX * SIN ( (RLI - 1.0) * ARG - PHI) 

ASIN(I) =-OMEG2 * XSIN(I) 



400 



CONTINUE 

COMPUTE THE EXPECTED VELOCITY AT 



THE MIDPOINT OF A DISPLACEMENT 



INTERVAL 

AMPV - XMAX * OMEGA 

ARG2 = ARG * 0.5 

VCOS(l) = AMPV * COS (ARG2 - PHI) 

DO 500 

RLI = I 
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VCOS(I) = AMPV * COS ( (RLI * ARG) - ARG2 - PHI) 

500 CONTINUE 

WRITE (3,16) OMEGA, XMAX 

16 FORMAT (' THE ANGULAR FREQUENCY AND AMPLITUDE HAVE VALUES' 2F10.2/) 
WRITE (3,18) 

18 FORMAT (' WHAT UNITS ARE ASSIGNED TO THE QUANTITIES LISTED ABOVE'/) 
WRITE (3,19) PHI 

19 FORMAT (' THE PHASE ANGLE HAS THE VALUE IN RADIAN MEASURE OF MINUS 
X • F10 .4 / ) 

WRITE (3,20) 

20 FORMAT ( ' POSITION RSIN (WT) VELOCITY RWCOS (WT) 

X ACCELERATION - RWWSIN (WT) ' ) 

WRITE (3,21) 

21 FORMAT ( ' CM CM CM/SEC CM/SEC 

X CM/SEC/SEC CM/SEC/SEC ' / ) 

WRITE ( 3,22) X(l), XSIN(l), V(l), VCOS (1) 

22 FORMAT (F10 . 2 , 6X, F9.2/ 31X, F9.1 , 6X, F9.1) 

WRITE ( 3,24) (X(I) , XSIN(I) , A(I) , ASIN(I) , V(I), VCOS (I) , 

XI = 2, Ji) 

24 FORMAT ( F10.2, 6X, F9.2, 36X, F9.0, 6X, F9.0/ 31X, F9.1,6X,F9.1) 
WRITE (3,26) X (J) ,XSIN (J) 

26 FORMAT (F10 .2, 6X,F9 .2) 

50 CALL EXIT 
END 
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APPENDIX 



Student Manual for Pendulum Experiment (Simple Harmonic Motion) 

The apparatus used in this experiment [Figure 1] is a pendu- 
lum made of a pointed metal bob suspended from a long wire. A 
strip of sensitized paper is placed parallel to the trajectory of 
the bob on a curved metal track. At regular intervals of time a 
large potential difference is applied between the bob and tne 
metal track causing a spark discharge. After one period of motion 
of the pendulum, the paper attached to the metal track contains a 
record of the position of the bob as a function of time. A typi- 
cal paper record is shown below the apparatus sketch [Figure 2] . 

We must extract from the tape not only position vs. time in- 
formation, but also the velocity and acceleration of the bob at 
equal intervals of time. This is a data reduction problem in which 
a large number of similar calculations must be performed. It is, 
therefore, an ideal problem to be handled by a computer. 

We assume that a computer is still a "black box" to most of 
you, i.e., that you have had little or no experience working with 
one. It is then appropriate to give a brief description of the 
function of a computer.* 

Imagine that you have employed a secretary who is efficient 
but not very intelligent. We can easily catalog the skills that 
she probably possesses. She can: 

1. perform simple mathematical operations — addition, 
subtraction, multiplication, and division. 

2. file information in definite locations and remember 
what these locations are. 

3. folJ.ow carefully detailed instructions step-by-step. 

4. compare two numbers and decide which is the greater. 

5. type information as it accumulates. 

These are not extreme standards to require in hiring a secretary. 
However, they represent the five functions performed by a computer. 
Why then are computers considered, for better or for worse, such 
an important part of our technological civilization? The answer 
is simple — speed. Because a computer can perform a calculation so 



*See, for example, Fortran for Physics , A. M, Bork, Addison 
Wesley Publishing Co., Inc., Reading./ Massachusetts (1967). 
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much more rapidly than any other calculating device, it can com- 
plete any problem in a brief period of time that a human could do 
only in an unreasonable amount of time. This change of time scale 
permits us to contemplate operations with a computer which would 
be in the realm of fantasy without this tool. 

Muhlenberg installed an IBM 1130 computer system ir< the fall 
of 1968. Its justification in an educational institution depends 
on its usefulness as a computational tool, as a bookkeeping device, 
and as a teaching machine in Skinner's sense. 

We wish to use the computer as a computational tool. To at- 
tack the problem of communicating with the computer, we must first 
decide what calculations we are to perform. Since we must make a 
graph of position, velocity, and acceleration and compare these 
with the theoretical values: 

X = R sin(ut - <p) 

V = v m cos (cat - <p) 

A - — sin ( cat — <p ) 

we see that we must calculate velocity, acceleration, and the quan- 
tities R, ca, and +. 

The program, i.e., the set of instructions which directs the 
computer to perform the calculations we need, is written in a lan- 
guage called FORTRAN (for FORmula TRANslation) . The language was 
designed to conform closely to the actual structure of algebraic 
calculations. A copy of the program is attached to this supplement. 
This program is stored in the memory of the machine; to run the 
program you need only enter data (the positions of the bob) in the 
form of punched cards. 

Perhaps the simplest way to understand both the computation 
procedure and the program is to look at the program and follow its 
statements through the calculation. Note that all statements pre- 
fixed by the letter C in the far left hand column are comments 
inserted into the program to clarify steps for the reader. They're 
ignored by the machine. We shall use them as "chapter headings.” 



C CALCULATE THE ANGULAR VELOCITY OF THE OSCILLATION 

The angular velocity u is equal to 2ir/T where T is the 
period of oscillation. We assume the positions of the bob (mea- 
sured from the midpoint of the motion) are entered into the machine 
consecutively beginning with the last negative value before the bob 
swings positive. We then find the time for which the pendulum bob 
had zero displacement. Subsequent positions are searched to find 
the next time the displacement was zero. The difference between 
these two times is one half the period; the angular velocity follows 
directly. 
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C CALCULATE THE PHASE ANGLE OF THE MOTION 

The displacement X is given by the expression R sin(wt-4>). 
When X is zero, the phase angle $ = «t. The statement after the 
comment referred to above performs this calculation. 



C CALCULATE THE APPROXIMATE AMPLITUDE OF THE MOTION 

This step illustrates quite nicely two programming techniques, 
use of the IF statement and use of the DO loop. 

The position data are denoted X(I) where I goes from one 
to the total number of data points. The technique we use to find 
the amplitude is simply to look at the absolute value of each X(I) 
and determine if it is larger than some other value of X(I) that 
we have assigned the value of XMAX. If it is larger we set the 
present value of X(I) equal to XMAX. XMAX is originally set 
equal to X(l) . 

Using the DO statement we let X(I) go from X(2) to X(J) 

(J is the total number of data points) by steps of one. Each 
value of X (I) is tested by the statement: 

IF (ABS (X (I) ) - XMAX) 100,100,90 
This means: 

1. Take the absolute value of X(I) . 

2. Subtract from this XMAX. 

3. If the result is 

a. negative — go to statement 100. 

b. zero — go to statement 100. 

c. positive — go to statement 90. 

Statement 90 changes XMAX from any previous value it had to a 
new value, viz. the absolute value of the present X(I) . The pro- 
gram then moves to statement 100 and then back to the beginning 
of the DO loop to increment I by one until I = J. 



C CALCULATE THE AVERAGE VELOCITY BETWEEN X(I) AND X(I+1) 

You have made this calculation in two previous experiments. 

V(I') = ( X(I+1) - X(I) )/ DELT 

where DELT is the time interval between two successive points. 

We identify this velocity as the average velocity evaluated at the 
midpoint of the interval between X(I) and X(I+1) . For this rea- 
son the velocity and position appear in the print-out spaced on 




alternate lines. 



C CALCULATE THE AVERAGE ACCELERATION BETWEEN V(I) AND V(I+1) 

This calculation is the same as that directly above 

A(I+I) = ( V(I+1) - V(I) )/ DELT 

represents the average acceleration between the two veloci- 
ties . 



The program concludes by calculating 

R sin (ut - 4>) 

Ru cos (ut - <f>) 

-Rw 2 s j * (cot - <f>) 

for the values of time associated with the computations above. 

Using the print— out, make three graphs: one of the displace - 

ment vs. time, the second of velocity vs. time, and the third of 
acceleration vs. time. Each should show experimental (or calcula- 
ted) and theoretical results . Why does the agreement between 
theory and experiment become steadily worse as we move successive- 
ly from displacement to velocity to acceleration? 



RELATIVISTIC TWO-BODY COLLISIONS 



Dale Winder 

Colorado State University 
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INTRODUCTION 



In the Physics Department at Colorado State University we 
have consciously endeavored to stimulate students to use the com- 
puter in their course work. In the introductory physics courses 
with calculus, we have fostered the generation of programs written 
by students and pertaining to the subject matter of the course on 
a voluntary basis* In the graduate- level Electromagnetic Theory 
course, problems were assigned to each student. The more experi- 
enced were given realistic problems, while novices were required 
to create function-generators. The program described here, ATOM, 
was written by a student in the third and last quarter of a course 
in general physics. It solves two typos of problems; 1) the rel- 
ativistic equations for the collision of two particles and 2) the 
transformation of coordinates between center- of -mass and labora- 
tory frame. It is unnecessarily complicated by including three 
dimensions and the possibility of conversion between mass and ki- 
netic energy. 

All participants in the introductory course were treated as 
individuals and asked to create a program on their own (with con- 
sultation and help, of course; . No manual was used, although a 
few introductory FORTRAN IV sessions were held for those who had 
no previous experience; in two one-hour sessions, I felt I could 
get enough across for the student to write simple programs involv- 
ing DO-loops, functions, data and printing. The chief virtue of 
this method was that the novice students could choose a subject 
which interested them, analyze it, program an extensive (for them) 
calculation and make it work. The next step, if we carry on with 
this program, would be to have them tackle a nontrivial calcula- 
tion, one that could not be computed in a finite time with a slide 
rule. The point in this effort is not in the quality of programs 
created, but in the involvement of the student in the use of a 
new, exciting, highly visible tool; one which many realize will 
be as common to them as a slide rule. There is an excitement and 
enthusiasm in solving problems of the relative difficulty and 
length represented by some of these problems. 
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STUDENT MANUAL 



The program, ATOM, described here solves the conservation 
equations for the collision of two sub-atomic particles. The input 
data consists of the mass, velocity and direction of motion of two 
particles before their collision and of one particle after the col- 
lision. The program then computes the mass, velocity and direc- 
tion of travel of the fourth particle. The total initial momentum 
of the two particles is then used to compute the transformation 
of variables to a frame of reference ("center-of -momentum" frame) 
where the total momentum of the system is zero. 

The fundamental physical principle involved is the conserva- 
tion of the relativistic momentum vector and the relativistic 
energy in the collision.^ For a single particle of rest mass mg/ 
speed v and velocity v these quantities are defined by the 
fundamental relativistic relations : 

m = m 0 (1 - v 2 /c 2 )“^ (relativistic mass) 
p = mv (relativistic momentum) 

E = me 2 (relativistic energy) 

which imply 

E 2 = p 2 c 2 + m^c 4 

relating energy, momentum and rest mass (c = speed of light) . 

For a two-particle system the conservation equations are 
thus expressed as 

Pi + P2 = P3 + P4 

Ei + E2 = E3 + E4 

assuming that two particles, 3 and 4, recoil from the collision, 
although they may not be of the same mass, individually or col- 
lectively, as the original colliding particles, 1 and 2. This is 
illustrated in Figure 1. 

For ease of programming, it is necessary to "scale" the 
variables appropriately; we shall do this by means of the follow- 
ing substitutions: 



BEFORE 



Z 

i u 
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/. f = pc = MV/y 
E = M/y and V = ?/E 

which also imply 0 

E 2 = P 2 + M 2 

The input data to ATOM are eight values Xj M l7 M 2 , M 3 , Vi, 
Vo, Vo, * t, and 0i and the index 1=1 for laboratory frame 
of reference, or 1=2 for center-of -momentum frame of reference. 
It is assumed that the direction of velocities and V determine 
the x-axis, and, although the 2 -axis could have been determined by 
requiring that it be in the plane of the x-axis and the recoil 
velocities, the program was purposely written in a general tnree- 
dimensional description. The direction angles, <j> and 0 were 
chosen as shown in Figure 1, and are assumed given m degrees. 



Every quantity to be output by the program is stored as a 
two-dimensional array of form X{I,N) where I specifies the 
frame of reference. The program starts by assigning storage space 
to MASS, speed (VEL) , magnitude of momentum (MOMENT) , JENERGY, * 
(PHI), 0 (THETA), vector components of v 3 (VEL^a) / Y k (VELi*^) / 
p 3r £ 4/ etc. In some arrays an extra space is set aside as a con- 
venient "working space." 



After specification of numerous FORMAT statements the program 
stores input data in the appropriate places (statement #41) , and 
proceeds (# 34 ) to compute momenta and energies for particles ±, 

3, and 4. After computing $i and 0i in radians, along with 
the appropriate direction cosines, the components of P 3 are cai 
culated (#707), and ? 4 is found fr °m the TOnscryation 
turn. $4 is computed from ? 4 (#23), followed by V 4 (#24+1) a 

finallv, the rest mass is found (#24+3) . The computation, for th 
first frame of reference, is completed with the calculation of the 
direction angles * 2 , 6 2 of the fourth particle, in degrees. 

Following this, the program transforms to the alternate frame 
of reference in ##702-800, essentially repeating the same computa- 
tions as in the preceding paragraph, after which the stored values 
are output (##36-40) in their appropriate formats. ^ listing ot 
ATOM and sample output is shown on the following pages. 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

17 

18 

19 

20 

50 

51 

52 

53 

54 
56 
58 



57 



505 

861 

862 

860 

41 

21 



34 



22 



PROGRAM ATOM 

REAL MASS (2,5) ,VEL( 2,4) , MOMENT (2,5) , ENERGY (2,5) , PHI (2.2) 

REAL THETA (2, 2) ,P(2*6),V(2,6),X(8),T(3) ,GAM(4) , RUTH (2) 

FORMAT (12,3 (F9 . 3) , 3 (F9 . 6) ,2 (F9 . 3) ) 

FORMAT (1H0,* INPUT DATA IS FROM LABORATORY FRAME OF REFERENCE*) 
FORMAT (1H0,* INPUT DATA IS FROM CENTER OF MOMENTUM FRAME OF REFERE 
XNCE*) 

FORMAT (1H0, 30X,*COLLISION AS SEEN FROM LABORATORY FRAME*, 13X,*COL 
XLISION AS SEEN FROM CENTER OF MOMENTUM FRAME*) 

FORMAT (1H0, 30X, 4 (*PARTICLE *,I1,2X, ) ,3X,4 (*PARTICLE *,I2,2X)) 
FORMAT (1H , 10X, *MASS (MEV) *, 7X,4 (F9 , 3, 3X) ,*XX*, 2X, 4 (F9 . 3, 3X) ) 

FORMAT (1H , 10X, * VELOCITY (X 1/C) *, 4X, 4 (F9 .3, 3X) , *XX*,2X, 4 (F9 . 3, 3X 
X) ) 

FORMAT (1H , 10X, *ENERGY (MEV) *, 7X, 4 (F9 . 3, 3X) ,*XX*, 2X, 4 (F9 . 3, 3X) ) 
FORMAT (1H , 10X, *MOMENTUM (MEV/C) *, 4X, 4 (F9 . 3, 3X) ,*XX*,2X, 4 (F9 . 3, 3X 
X) ) 

FORMAT ( 1H0 , 1GX , *MOMENTUM* ) 

FORMAT (1H , 14X, *COMPONENTS* ) 

FORMAT (1H / 19X, *X*, 10X, 4 (F9 .3, 3X) ,*XX*, 2X, 4 (F9 . 3, 3X) ) 

FORMAT (1H / 19X, *Y*, 10X, 4 (F9 . 3, 3X) ,*XX*,2X, 4 (F9 .3, 3X) ) 

FORMAT (1H / 19X,*Z*, 10X, 4 (F9 .3, 3X) ,*XX*, 2X, 4 (F9 . 3, 3X) ) 

FORMAT ( 1H0 , 10X, *VELOCITY* ) 

FORMAT (1H ,19X,*X*,10X,4 (F9.6,3X) ,*XX*,2X,4 (F9 .6, 3X) ) 

FORMAT (1H / 19X, *Y*, 10X, 4 (F9 .6, 3X) ,*XX*,2X,4 (F9 .6, 3X) ) 

FORMAT (1H , 19X, *Z*,10X, 4 (F9.6,3X) ,*XX*,2X,4 (F9.6, 3X) ) 

FORMAT (1H0 , 10X, *DIRECTION OF*) 

FORMAT (1H , 10X , *TRAVEL* ) 

FORMAT (1H ,13X,*PHI (DEGREES) *, 4X, 4 (F9 . 3, 3X) ,*XX*, 2X, 4 (F9 . 3, 3X) ) 
FORMAT (1H , 14X, *THETA (DEGREES) *,4 (F9 .3, 3X) ,*XX*,2X,4 (F9 . 3, 3X) ) 
FORMAT (1H0 , 10X, *TOTAL MOMENTUM*, 2 4X,F 9 . 3, 21X,*XX*,20X,F9 . 3) 

FORMAT (1H , 10X/*TOTAL ENERGY *,24X,F9.3,21X,*XX*,20X,F9 .3) 

FORMAT (1H0) 

FORMAT (1H0 ,20X, ‘VELOCITY OF CENTER OF MOMENTUM AS SEEN FROM LABOR 
XATORY = *,F18. 7,* TIMES C .*) 

ZZ = 57.2957795 
A1 = 0.0 
A2 = 90.0 
A3 = 180.0 
A4 == 0.0 

READ (5, 1) (I, (X (J) , J=l, 8) ) 

GO TO (39, 40) , I 

IF (X (4) ) 861, 862, 862 
A4 = 180.0 
IF (X (5) ) 41, 860, 860 
A3 = 0.0 
DO 21 J = 1,3 
MASS (I, J) = X(J) 

VEL(I, J) = X (J+3) 

PHI (1,1) = X (7) 

THETA ( I, 1) = X (8) 

DO 22 J = 1,3 

GAM(J) = SQRT(1.0 - VEL(I,J)**2) 

MOMENT (I, J) = (MASS (I, J) * VEL(I,J)) / GAM(J) 

ENERGY (I, J) = MASS (I, J) / GAM ( J ) 




MOMENT (1,5) = MOMENT (I, 1) + MOMENT(I,2) 

ENERGY (1,5) = ENERGY (I, 1) + ENERGY(I,2) 

ENERGY (I, 4) = ENERGY (1,5) - ENERGY (I, 3) 

PHIRAD = PHI (1,1) / ZZ 
THETAR = THETA (1,1) / ZZ 

T (1) = SIN (PHIRAD) * COS (THETAR) 

T (2) = SIN (PHIRAD) * SIN (THETAR) 

T (3) = COS (PHIRAD) 

DO 707 J = 1,3 
V(I,J) = VEL (I, 3) *T (J) 

707 P(I,J) = MOMENT (1,3) *T ( J ) 

P (1,4) = MOMENT (1, 5) - P(I,1) 

P (1,5) = - P (1,2) 

P (1,6) = - P (I, 3) 

M = 4 

DO 23 J = 4,6 

23 V(I,J) = P(I,J) / ENERGY (X,M) 

SUMl =0.0 

SUM2 =0.0 

DO 24 J = 4,6 

SUMl = SUMl + V (I, J) **2 

24 SUM2 = SUM2 + P(I,J)**2 
VEL (I, 4) = , SQRT (SUMl) 

MOMENT (I, 4) = SQRT (SUM2) 

MASS (1,4) =SQRT (ENERGY (1,4) **2-MOMENT (1 , 4 ) * *2 ) 

X = ACOS (V (I, 6)/VEL (I, 4) ) 

PHI (1,2) = X * ZZ 

XX = ATAN ( V ( 1 , 5 ) /V (1,4)) 

THETA (I, 2) = XX * ZZ 
IF ( V ( I , 5 ) ) 701,702,702 

701 THETA (1,2) = 360.0 - THETA (1,2) 

702 IF (I .EQ.2) GO TO 26 
II = 2 

DO 704 J=l, 2 

704 RUTH (J)= MASS (I, J) /GAM (J) 

VCM= (RUTH (1) *VEL (I, 1)+RUTH (2) *VEL (1,2) )/ (RUTH (1)+RUTH (2) ) 
GO TO 27 

26 II = 1 

VCM= -VEL (1,2) 

27 DO 28 J=l, 4 

28 MASS (II, J) = MASS (I, J) 

DO 29 J=l, 2 

29 VEL (II, J)= (VEL (I, J) -VCM) / (1. 0-VCM*VEL (I, J) ) 

GAMMA=SQRT ( 1 . 0-VCM**2 ) 

DO 31 J=l, 4, 3 
V7=1.0- VCM*V (I, J) 

V(II,J) = (V (I - J) —VCM) /V7 

K1=J+1 

K2=J+2 

DO 31 L= K1,K2 

31 V (II, L) = (V(I,L) * GAMMA) / V7 
SUMl =0.0 
SUM2 =0.0 
DO 32 J=l, 3 
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32 



33 



25 



727 



800 

36 



SUM! = SUM1 + V(II / J)**2 
Kl = J+3 

SIJM2 = SUM2 + V (II, Kl) **2 
VEL (II, 3) = SQRT(SUMl) 

VEL (II, 4) = SQRT (SUM2) 

DO 33 J=l, 4 

GAM(J) = SQRT (1.0 - VEL (II, J) **2) 

ENERGY (II, J) = MASS (II, J) / GAM(J) 

MOMENT (II, J) = (MASS (II, J) * VEL (II, J) ) / 

MOMENT (II, 5) = MOMENT (11,1) + MOMENT (11,2) 

ENERGY (II, 5) = ENERGY (11,1) + ENERGY (11,2) 

M=3 

DO 25 J=l, 6 

V7 =SQRT (X. 0— V (II, J)**2) 

IF (J.EQ. 4) M=4 

P(II,J) = (MASS (II, M) * V(II,J)) / V7 

L7 = 3 

K = 3 

N = 2 

N4 = 1 

DO 36 J=l,2 

IF (J.EQ. 1) GO TO 727 

L7 - 4 



K = 6 
N = 5 
N4 = 4 

X = ACOS (V (II,K) / VEL (II, L7) ) 
PHI (II, J) = X * ZZ 
XX = ATAN (V (II, N) / V (II, N4) ) 
THETA ( II, J) = XX* ZZ 
IF (V (II,N) ) 800,36,36 
THETA ( II, J) = 360.0 - THETA (II, J) 
CONTINUE 
WRITE (6, 4) 

WRITE (6,5) 

WRITE (6, 6) 

WRITE (6, 7) 

WRITE (6, 8) 

WRITE (6, 9) 

WRITE (6, 10) 

WRITE (6, 11) 

WRITE (6, 12) 

WRITE (6, 13) 

WRITE (6, 14) 

WRITE (6, 15) 

WRITE (6,11) 

WRITE (6, 17) 

WRITE (6 18) 

WRITE (6, 19) 

WRITE (6,20) 

WRITE (6,50) 

WRITE (6, 51) 

WRITE (6,52) 



((I, 1=1, 4), (J,J=1,4)) 

( (MASS (I, J) , J=l, 4) ,1=1,2) 

( (VEL (I, J) , J=l, 4) ,1=1,2) 

( (ENERGY (I, J) , J=l,4) ,1=1,2) 
( (MOMENT (I, J) , J=l,4) ,1=1,2) 



( ((MOMENT (J, I), 1=1, 2) ,P(J,1) ,P(J,4) ) ,J=1,2 
( (A1,A1,P(J,2) ,P (J, 5) , J=l,2) ) 

( (A1,A1,P(J,3) ,P (J, 6) , J=l,2) ) 



( ((VEL(J,I) ,1=1,2) ,V(J,1),V(J,4)) , J=l,2) 
( (A1,A1,V (J ,2) ,V(J,5) , J— 1,2) ) 

( (A1,A1,V(J,3) ,V(J,6) , J=l,2) ) 



( (A4,A3, PHI ( J , 1 ) , PHI ( J , 2 ) ,J=1,2)) 

( (A2,A2, THETA ( J, 1) ,THETA(J,2) , J=l,2) ) 




GAM(J) 




WRITE (6,53) ( (MOMENT ( J , 5 ) ,J=1,2)) 

WRITE (6, 54) ( (ENERGY (J, 5) , J=l,2) ) 

WRITE (6, 58) ( VCM ) 

DO 55 J=l, 5 

55 WRITE (6 ,56) 

GO TO 57 

39 WRITE (6,2) 

GO TO 505 

40 WRITE (6, 3) 

GO TO 505 
END 
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TEACHER* S GUIDE 



Available for student use were a CDC 6400 and an IBM 1401, 
neither one in a time-sharing mode; the language was FORTRAN. Job 
cards were issued upon request by the instructor at the beginning 
of the quarter, each good for one ten-second, ten-page run on the 
6400 at a cost of $1. Students prepared, submitted and debugged 
programs at their convenience. A "pep talk" on the importance and 
relevance of computer calculations in science and engineering and 
a selection of possible problems helps to involve the students and 
to assure participation. 

This program, ATOM, has run well in numerous tests, although 
it is not safeguarded against zeros in denominators or imaginary 
roots which could occur, because input is so free. In most cases, 
we have restricted ourselves to elastic collisions in a plane (0i, 
= 0, M 3 = M x or M 3 = M 2 ).* 



* Editor*s Note : It might also prove of interest to compare 

output from this program with that from the nonrelativistic pro- 
gram MOMEN in D. T. Grimsrud*s "Two Experiments" in Computer - Based 
Physics : An Anthology , R. Blum, et al, published by the Commis- 

sion on College Physics (1969) . Since conservation computations 
are performed in terms of momentum and energy, it should be possi- 
ble to modify ATOM to include Compton scattering as a special case 
for which, say, V 1 =1=V 3 . A logical branch in the case of in- 
put V* = 1 would then compute E^ = Mj = hv = hc/X and P^ = E^ 
for the photon, omitting such undefined computations as #34-#22, 
etc. 
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l ERIC 



INTRODUCTION 



A discussion of high energy particle physics is an integral 
part of most modern physics courses. Laboratory exercises using 
bubble chamber photographs are commonly employed to amplify the 
experimental side of these studies while questions of conservation 
laws and allowed reactions are covered in class and problem as- 
signments. It is often very difficult for a student to capture 
the true flavor of this field from classroom treatments alone. 

Can a student be led to discover, or at least to test, the con- 
servation laws as an active research participant? Is there a way 
for him to build the particle spectrum as a result of his own 
ideas and imagination? 

The computer-based laboratory exercise presented here is a 
first attempt to answer these questions. In this exercise the 
student is asked to imagine himself as a theoretical physicist who 
suggests experiments to a large accelerator laboratory. The gen- 
eral idea is to allow the student to predict high energy particle 
reactions and to test his predictions by " experiments " carried out 
in the computer- simulated accelerator laboratory. 

This laboratory exercise has been used during one term in a 
pre-calculus course for freshmen, mostly prospective physics and 
mathematics majors, called Introductory Modern Physics.* Since 
the classroom discussion of the concepts involved in the laborato- 
ry exercise did not come until the last week in the term, these 
students were really given an opportunity to experience life as it 
confronts the working physicist. Some took the exercise as a 
challenge and only went to the literature after exhausting their 
ingenuity and understanding. Others began by reading about the 
field and then proceeded to test and verify their understanding. 
Either approach seemed rewarding and succeeded in getting the stu- 
dents to study the field with real enthusiasm. The only computer- 
related skill required of the students for this exercise was card 
punching c As used in this course the computer experiments by the 
students were carried out during a six week period, so a fast turn- 
around time was not essential. The computer used by the class was 
an IBM 1130 operated in a hands-on, open shop where the students 
could either submit jobs for batch processing or operate the com- 
puter themselves. 

A revised version of the "Student Manual" that was used by 
38 students in the fall of 1968 is presented in the next section 
which will both illustrate the flavor of the exercise and give the 



* Interested parties may write the author to request further 
information about this unique course. 
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reader an introduction from the user's point of view. The "Teach- 
ers Guide" which follows it contains more detailed information 
about the simulation program and its use by students and concludes 
with the data storage techniques used with the IBM 1130 system, a 
listing of the particles and reactions contained in the simulation 
and complete program listings. 



STUDENT MANUAL 



Introduction 



One of the most interesting fields of contemporary physics 
is particle physics, sometimes called "high energy" physics, be- 
cause the paftieles involved in the experiments have very high 
kinetic energies. The reactions take place in the GeV (10 9 eV) 
region and above, whereas chemical reactions involve energies on 
t ; e order of one eV or so. 

Now the question comes to mind, "How can those of us at a 
liberal arts college, without access to a high energy accelerator, 
get some firsthand knowledge of this field?" In our course on 
modern physics we attempt to attack this problem by including two 
experiments related to this field. One of these has become an 
old standard in the course, the analysis of some actual bubble 
chamber photographs.* A rather extensive write-up is available 
for that experiment and you should read that write-up before pro- 
ceeding with this new experiment. It is not necessary for you to 
have actually worked on the bubble chamber photographs before 
attempting this new project; but, it might make things clearer if 
you have. 

The new experiment, the subject of this write-up, requires 
that you think of yourself as a theoretical "high energy" physi- 
cist. You are working at one of the large accelerator laborator- 
ies where you help plan new experimental projects and try to make 
sense out of the experimental results. This allows you to ponder 
the results of experiments, to ask the questions that lead to new 
experiments and, hopefully, to discover the answers to your ques- 
tions . 

You might well ask, "What facilities do I have at my command?" 
The answer to that is, "An entire accelerator laboratory — acceler- 
ator, hydrogen bubble chamber, technical staff, and a crew to scan 
the bubble chamber pictures." It will be up to you to decide on 
the type of particle you would like to fire into the bubble cham- 
ber. You will also choose the kinetic energy of this incident 
particle. And, finally, you get to instruct the scanners to search 
for up to three types of resulting particles by simply specifying 
the masses of these "predicted" particles. This is all accom- 
plished by letting the computer simulate the entire operation of 
the accelerator laboratory. 

By studying the types of reactions that can, and do, take 
place at very high energies, a list of "particles" was selected as 
being representative of those actually known to physicists. The 
initial configuration of the simulation allows only protons to be 



*See Welch Catalog, No. 2171 — Bubble Chamber Photo Analysis 
Apparatus . 
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the target particles (a liquid hydrogen bubble chamber, remember?) . 
Initially, the incident projectile particles must be ir “ mesons, 
and as work progresses on the simulation other projectile parti- 
cles will become available; but you have to discover them first 
using the tt ~ mesons. The work involved here is not just getting 
the data — the laboratory supporting personnel do that for you — 
but rather in trying to make some sense out of it and planning new 
experiments . 

If you want to know more about the workings of an actual ac- 
celerator and how a it" + p reaction comes about, you should read 
the write-up for the bubble chamber photo analysis experiment. If 
you would like to know more about the particles you will "discov- 
er" during the course of your experiments, you should read Chapter 
10 in Beiser*s paperback Modern Physics ; An Introductory Survey , 
Addison- Wes ley, Reading, Mass. (1968) , and from there head for the 
library to the latest books on nuclear physics and pair tide phys- 
ics, and to Scientific American articles. A list of possible 
references may be found in the book by M. S. Livingston, Particle 
Physics , McGraw-Hill, New York, N.Y. (1968), pp. 222-224. If you 
want suggestions for further study, please ask your instructor. 

As a theoretical physicist, what parameters do you think are 
necessary to specify a particle, that is, to differentiate it from 
another particle? Suppose you are told there exist two particles 
with rest mass of about 939 MeV/c 2 , (c = speed of light) one with 
positive charge and a second with no charge. At the present time, 
physicists refer to these two particles as two "charge states of 
the nucleon." You already know these particles by the names pro- 
ton and neutron. In addition we know of the anti-nucleon states, 
the anti-proton and anti -neutron. We thus use charge as one means 
for differentiating particles. As you read about the elementary 
particles, you will discover many other ways for differentiating 
between them, but keep in mind that particles of equal mass but 
different charge may very well be simply different charge states 
of the same "fundamental particle" like the "nucleon." 

As you read Chapter 10 in Beiser’s book, you will become 
aware that the lifetime of the particles is an important parameter. 
Reactions that are "allowed" to take place by a strong interaction 
will happen in the short time it takes for two particles to pass 
one another while moving at a relative speed of approximately c. 

If we take the diameter of a nucleon (proton or neutron) to be 
about 10“ 13 cm (one "fermi") this time comes out to be about lO -24 
seconds. However, Mother Nature hasn’t "allowed" everything to 
happen that fast. In nuclear physics some reactions are called 
"weak" interactions and they are found to take place about 10 
times slower than the strong interactions. The time associated 
with weak interactions is thus about 10“ 10 seconds. You will find 
both strong and weak interactions going on as indicated by the 
lifetimes . 




In the procedure section you will be guided through a trial 
run. It is suggested that you have the computer do this run for 
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you as the first step in using the "laboratory.” The laboratory 
will be at your disposal during the remainder of the term so you 
may feel free to do experiments any time. It is intended that 
this exercise occupy at least one week's worth of laboratory time. 
Be sure to read the note at the end about "publishing" your re- 
sults . 



Experimental Procedure 

For each experimental run you must supply the accelerator 
staff with the identity of your projectile particle (a tT meson 
initially) by specifying its mass and charge (137 MeV/c 2 and -e 
for ir”) . Then you must specify the laboratory kinetic energy of 
this incident particle in MeV; 0-2500 MeV will cover the entire 
range of simulated reactions. Only reactions with threshold en- 
ergies within 10% of your specified energy will be detected. 

Finally, once you decide on a projectile particle, you must 
tell the picture scanners what you are looking for. That is, you 
must "predict" the masses of up to three particles you expect to 
have produced. The overwhelming number of pictures obtained from 
an experimental run prohibits a complete analysis, so the photo 
scanning process is programmed to reject any reactions that do not 
result in one of the three particle masses specified by you. The 
scanners can determine masses within ±5%. 

As a theoretician you are interested in the results of the 
experiments. When the scanning crew finds a particle with a mass 
specified by you, they will record the data and give you a run- 
down of the significant reactions that occur during the experiment. 
For example, suppose the it” + p collision produces two parti- 
cles: (1) mass Mi and charge Qi , and (2) mass M 2 and charge 

Q 2 . You will be told these masses and charges. Suppose further 
that particle (2) decayed after 10" 8 seconds into two other par- 
ticles that also appeared on the photograph. You will also be 
told the masses and charges for these decay products and the life- 
time of the parent particle (2) . Since these aie’ chreshold reac- 
tions -the reaction products are assumed stationary in center-of- 
mass coordinates . 

In order to get down to earth let us consider the ir + p 
reaction as a trial run. Do you suppose there exists a particle 
with half the nucleon mass? To be specific, we want to find out 
whether the reaction tt” + p N + M(500) exists. Here we have 
written N to stand for a nucleon (regardless of its charge state; 
i.e., either a neutron or proton) and M(500) stands for the un- 
discovered particle with rest mass around 500 MeV/c 2 . The in- 
crease in mass-energy involved in producing the 500 MeV/c 2 parti- 
cle is the difference between the new mass and the pion mass, 
500-137 = 363 MeV/c 2 . But we also know that some of the pion kine- 
tic energy is effectively "lost" as center -of -mass energy and is 
not available for appearance as rest mass. Using the conserva- 
tion of momentum and total energy we can calculate the threshold 
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kinetic energy required for the pion to be able to produce 

this new particle. The result of this derivation is 

Eth (threshold kinetic energy) = QM/2m 2 (Eq. 1) 

where: Q = (initial total rest mass - final total rest mass)c 2 

M = (total final plus initial rest masses) 

m 2 = (rest mass of target particle) 

To be specific, Q = (939 + 137 - 939 - 500) MeV = 363 MeV 

M = (939 + 500 + 939 + 137) MeV/c 2 = 2515 MeV/c 2 

m 2 = 939 MeV/c 2 

Carry out the arithmetic and you predict a threshold kinetic ener- 
gy of about 485 MeV as required for the incident pion. Let*s pick 

500 MeV as a round number recalling that our scanners give us a 

5% tolerance on mass determinations. 

The FORTRAN READ statement for this program is the following: 

READ (2,201) Ml, IQ1, KE, MOUT(l), MOUT(2), MOUT(3) 

201 FORMAT (615) 



For input data we specify, as integers, the incident particle mass. 
Ml, its charge in units of e, IQ1, its kinetic energy, and than 
up to three mass values M0UT( ) we want the scanners to look for. 
If we want the scanners to look for the new 500 MeV particle, a 
nucleon and a pion, the data card would take the following form: 



137 -1 500 500 939 137 

In addition we need some XEQ cards and cards to define some stor- 
age areas. The cards required to call the program and to specify 
our data are as follows: 



Columns : 



5 10 15 20 25 



30 35 40 



Card Is 
Card 2 : 
Card 3 : 
Card 4: 
Card 5 : 



// JOB 

// XEQ PHYHE 01 

♦FILES (10, PARTS) , (1,REAC1) , (2,REAC2) 
137 -1 500 500 939 137 

blank 



The computer will digest these five cards and print the following 
messages : 
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INITIAL EXPERIMENTAL CONDITIONS 

INCIDENT PARTICLE MASS= 137 MEV/C2 WITH CHARGE= - 
KINETIC ENERGY OF PARTICLE BEAM = 500 MEV 



I 

i 

I 



SCAN IS FOR MASSES 500, 939, AND 137 MEV/C2 

** PI MESONS ARE TOO NUMEROUS FOR A SCANNING REPORT 
THE SCANNERS ARE IGNORING YOUR REQUEST 



1 



SCAN HAS DETECTED NO PARTICLES WHICH SATISFY YOUR CONDITIONS. 



As you can read, the computer will verify that it understands your 
specified initial conditions and then repeats your request for the 
scanning masses. In this case the scanners will not look for pi 
mesons (MOUT(3)) because they are simply too numerous, and finally 
we are told that no reaction of the kind we sought has been found. 

Now here is the first opportunity to make an educated guess. 
Suppose you have a theory that some sort of particle exists which 
is really a collection (molecule?) of five bound pions. You might 
guess its mass to be 5x137=685 MeV/c 2 . If you use this mass to 
calculate the laboratory threshold kinetic energy (Eq. 1) you pre- 
dict about 790 MeV as the required energy for the incident pion. 

We can now do a new experiment simply by replacing Card 4 by a new 
card. 

Card 4: 137 -1 800 500 939 685 

where the incident pions have a kinetic energy of 800 MeV and the 
scan will be made for particles of mass 500, 939, and 685 MeV/c 2 . 
The program yields the output shown on the next page (Figure 1) . 

Surprise! There ijs a particle with about half the mass of a 
nucleon. But, the most surprising part seems to be that it is not 
produced as simply as we had supposed. We had to supply more en- 
ergy, apparently enough so that a new heavier particle could be 
produced to accompany our 496 MeV particle. Why didn’t the neu- 
tral 496 MeV particle appear in the first experiment when there 
was certainly enough energy? Is there a conservation law acting 
here that prevents the production of the neutral 496 MeV particle 
in conjunction with, say, a neutron? Are there also charged par- 
ticles with the masses 1115 MeV/c 2 and 496 MeV/c 2 ? We should note 
that no 685 MeV/c 2 particle was found. Does that mean it doesn’t 
exist or that we simply don't know how to produce it? 

If you look at the decays of these two new particles you find 
that the 1115 MeV particle decays in 2x10“ 10 seconds into either 
a proton and pi-minus meson or into a neutron and pi-zero meson. 

A weak interaction? What prevents a fast decay? Now look at the 
496 MeV particle's decay habits. Nothing but puzzles! It seems 
to be a Jekyll-and-Hyde particle. Is it one particle with two 
different decay paths, or two particles with the same mass and 
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INITIAL EXPERIMENTAL CONDITIONS 

INCIDENT PARTICLE MASS= 137 MEV/C2 WITH CHARGE= -1 

KINETIC ENERGY OF PARTICLE BEAM = 800 MEV 

SCAN IS FOR MASSES 500, 939, AND 685 MEV/C2 



SUCCESS. RESULTS FOLLOW. 



PARTICLE (1) 
PARTICLE (2) 



CG= 0 
CG= 0 



MASS= 

MASS= 



1115 

496 



DECAY TIME = 
DECAY TIME = 
DECAY TIME = 



PARTICLE 


(1) FOLLOW. 


9 




939 


0.20E-09 


CG = 


1 


MASS 


— 




CG = 


-1 


MASS 


zz 


137 


0.20E-09 


CG = 


0 


MASS 


= 


939 




CG = 


0 


MASS 


— 


137 


PARTICLE 


(2) FOLLOW 


• 






0.70E-10 


CG = 


1 


MASS 


— 


137 




CG = 


-1 


MASS 




137 


0.70E-10 


CG = 


0 


MASS 


= 


137 




CG = 


0 


MASS 


zz 


137 


0.40E-07 


CG = 


1 


MASS 


zz 


137 




CG = 


-1 


MASS 


zz 


137 




CG = 


0 


MASS 


zz 


137 


0.40E-07 


CG = 


0 


MASS 


zz 


137 




CG = 


0 


MASS 


zz 


137 




CG = 


0 


MASS 


zz 


137 



Figure 1 
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charae? It either decays into two pions with a lifetime of 0.7 
xlO-I° seconds or it decays into three pions with a longer life- 
time of 4xl0~ 8 seconds. What is going on here? 



Conclusion 



Now you are on your own. It is hoped that you will come to 
appreciate the methods, joys, and tribulations of the theoretical 
physicist and that you will learn about the strange world of high 
energy physics, while having fun doing it. You should not sit in 
the computer room grinding out lots of numbers; instead, try to 
think up experiments that will provide you with answers to your 
questions. Keep track of your hypotheses, experiments, and mis- 
cellaneous thoughts in your lab notebook. The reactions avail- 
able to you do not exhaust all of those possible in nature, but 
every effort has gone into providing a representative sample. 

If you are lucky enough to discover new particles, you should 
"publish" your results on the bulletin board. Your instructors 
have started a "publication sheet" with their discoveries; you can 
add your results to it, unless someone else publishes first. How- 
ever, verify ideas carefully before "publication," and beware the 
occupational disease of R.I.P. (Rushed Into Print). 

To communicate your guesses or questions to other class mem- 
bers, simply write a note and stick it on the bulletin board. It 
is hoped that we may develop a small community of researchers, all 
working and sharing results in their common goal of understanding 
nature . 



o 



TEACHER'S GUIDE 



Many simulations are designed for real time interaction be- 
tween the user and the simulator. This acceleration laboratory 
simulation is not. Rather, it is intended that a student will 
try experimental runs primarily as tests for his ideas and will 
spend some period of time between these runs thinking about the 
results and planning new experiments. In order to encourage this 
and to discourage a shotgun approach the program presented here 
reads only one experimental specification per program execution, 
hence, it is readily usable with nearly all computer systems, 
time— shared or batch processed. In some cases there may be cost 
considerations, a desire for short turnaround times, or limita- 
tions on the number of computer runs allowed by any individual or 
class per day. In situations like these, being able to read more 
than one experimental specification per computer job could be 
easily achieved through a minor modification of the source pro- 
gram, such as relying upon an unsatisfied READ statement or a 
test for negative Ml to terminate the program. 

Our students were encouraged to ’publish" their discoveries 
by filling in a chart that was tacked on a bulletin board outside 
the lecture room. The chart included columns for particle charge, 
mas s, lifetime, decay paths, and production technique. Additional 
space on the bulletin board was reserved for theoretical "papers" 
which consisted of short notes (by both students and instructor) 
asking questions about the discoveries and pointing out special 
properties. The students who worked hardest on the experiments 
quickly discovered many of the easy particles and the other stu- 
dents required reassurance that not everything had been discover- 
ed. One pair of students discovered that the K“ was available 
for use as an additional incident particle before this fact was 
revealed to the class by the instructor. Once the K became 
available many new reactions were possible and the late starters 
were able to get into the act. 

It would be a mistake to expect students to arrive at a com- 
plete understanding of all the variations possible in this exer- 
cise without a fair amount of outside help. In addition to being 
encouraged to share their results and ideas with each other, our 
students were expected to use other references, the two primary 
ones being the book by Livingston, Particle Physics , McGraw-Hill, 
New York, N.Y. (1968) and the article by Chew, Gell-Man, and Ro- 
senfeld. Scientific A merican (February 1964) .* 



*See Editor's Note on following page. 



96 



An opportunity for classroom discussion is essential for 
collecting loose ends. Coming after the experimentation, such a 
discussion can draw examples from the students" own "published" 
results and can point out the limitations in the simulation 
program. Furthermore, having prior information about the baryon 
spectrum did not seem to dampen enthusiasm among our students. 
Indeed, a theoretical prediction that a certain particle cannot 
be produced by either the ir” + p or the K~ + p reaction can 
be of more value than the publication of an accidental discovery 
of seme new particle. 

This simulation exercise was developed for use in the lab- 
oratory portion of an introductory course in modern physics. 
However, there do seem to be other possible uses for this pro- 
gram. For example, if a time-shared terminal could be made 
available in the classroom, such a simulation program might be 
very useful in a discussion session. Another possibility might 
be a problem set where the students could use the simulation in 
response to leading questions or as a check on their answers. 



Editor" s Note: To obtain the formula Eth = QM/2m 2 given 

in the Student Manual, recall that for an assembly of particles 
the quantity 

(total energy) 2 - (total momentum) 2 c 2 

is invariant with respect to Lorentz transformations and that en- 
ergy is conserved in the collision. In the center-of-mass system 
the total momentum is, by definition, zero; and the energy before 
collision, E', equals E", the energy after collision. However, 
from the invariant, E' and E" are related to laboratory obser- 
vations as follows: 



; * 2 — - 4- Ttt ~ -}- P.-uU ^ 2 

! 4 *»* / * J * 



Pl 2 c 2 ; 



E" = M" 



where M" is the rest mass of the reaction products and pi the 
relativistic momentum of the projectile particle 

Pl 2 c 2 = (m x + Eth) 2 “ ™l 2 ~ 2m i Eth + E th 2 

Making the appropriate substitutions in the equation E* 2 = E" 2 
yields 

M" 2 - (m! + m 2 ) 2 = 2Ethm 2 



and factoring the left-hand side of the equation as the differ 
ence of two squares leads directly to the desired relationship. 

This method is discussed, with examples, in Appendix A of 
Introduction to Elementary Particles, by W. S. C. Williams, Aca- 
demic Press, New York, N.Y. (1961) . The result can also be de- 
rived by straightforward but tedious calculations from conserva- 
tion of energy and momentum. These can be found worked out in 
detail in The Atomic Nucleus , by R. D. Evans, McGraw-Hill, New 
York, N.Y. (1955) , Chapter 12. 




Many of the limitations contained in the present program are 
artificial and were included in order to reduce the development 
and programming time or to limit the student* s freedom. For ex- 
ample : 

1. The only particles available in the program are the 
baryons up to mass 1688 MeV/c 2 , the 7r and K mesons, 
and the muons. 

2. The particle manifold masses are combined. For ex- 
ample, both the proton and neutron appear as having a 
mass 939 MeV/c 2 and all the ir-mesons appear with mass 
137 MeV/c 2 . 

3. Only reactions having thresholds within 10% of the 
specified threshold energy are detected. This pre- 
vents picking a high energy and learning about all the 
reactions with lower thresholds in a single experiment. 

4. No possibility has been included for reactions pro- 
duced by secondary particles. Cnly 7r“ + p and K“ 

+ p reactions and the decays of their resulting parti- 
cles are included. 

5. No information on the relative frequency of the vari- 
ous reactions and decays is included. 

6. Electrons, gamma rays and neutrinos are ignored as 
being undetected. 



A good case can be made for eliminating many of these limitations 
in future versions of the simulation, particularly if it is to be 
used by advanced students. It would not be difficult to include 
more members of the meson family or to increase the energy range 
and number of available reactions. The decay schemes could be 
greatly improved if electrons and gamma rays were included in the 
particle list. 

The Program 

The basic operation of the program is as follows : 

1. Input data is read and a list (Table 1) of 33 particles 
is scanned to see if any of the given masses correspond to items 
on the list. Each item is indexed and the identifying indices of 
legitimate particle masses stored for later reference. 

2. The program then scans a list (Table 2) of 24 possible 
7r” + p or 32 possible K~ + p reactions, depending on Ml and 
the threshold energy to determine which, if any, could give rise 
to the anticipated particles. 

3. If an otherwise eligible reaction from Table 2 is found 
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not to contain the anticipated particles, then Table 1 is checked 
again to see if they might not be decay products of the reaction 
being examined. If successful, the program writes the result be- 
fore returning to the next eligible reaction, if any, in Table 2. 

As an illustration of how to read the particle table. Table 
1 on the following page, we shall consider particle 21 which has 
a mass of 1405 MeV/c 2 , charge 0, mean life 0.5x10“ 2 3 seconds. 

The 12 numbers following the lifetime are the decay paths record- 
ed in four groups of three particles. It follows that particle 
21 has three decay paths (zeros have no meaning) into the parti- 
cle pairs (9,-1), (10,2), and (11,1). These integers refer to 
other particle numbers and may be looked up in the particle table. 
Negative particle numbers indicate the anti -particle. From the 
last column of Table 1 we can identify the particle as a A 0 *? 
its three decay paths are (l+,ir“), (I°,ir°), and (l“,ir + ) . 

In Table 2 (page following Table 1) the reactions are listed 
in order of increasing threshold energy. Reaction products are 
j given by particle numbers as found in the particle table. For 

example, consider reactions 13 and 14 having a threshold energy of 
j 219 MeV. These produce particles (7,4,2) and (7,5,2), which may 

j be identified as being (n,Ki°,ir°) and (n,K 2 °,ir°) . 

This simulation program was developed for use with the IBM 
1130 version 2 disk monitor system. In this system, the above 
particle table and reaction lists were stored as DATA FILES on 
the system disk, using the auxiliary program listed in Figure 2. 

Some readers may not be familiar with the DEFINE FILE state- 
ment or with the Disk WRITE (NUM'J) XXX statement which appears as 
statement numbers 1 and 2 in this program. The DEFINE FILE 
statement gives information allowing the compiler to set up a 
data file on the disk. The WRITE ( * ) statement is an instruc- 
i tion to "write on disk." Further details may be found in Using 

the IBM 1130 , by A. Bork, Addison-Wesley , Reading, Mass. (1968). 

The simulation program itself was also stored on the disk 
•and could be called by the students using the short five- card 
deck listed in the Student Manual. This program was stored under 
the name PHYHE (Physics High Energy) , hence the presence of that 
name on the program execute card (// XEQ) . The third card of 
that deck is used to relate files numbered 10, 1, and 2 in the 
program READ( 1 ) statements to the names PARTS, REAC1, and REAC2 
that were given to the particle table and reaction lists when 
stored on the disk. The fourth card contains the input data, and 
the last blank card simply insures that the card reader will read 
the fourth card if no jobs follow this one. 

The flow chart on page 100 is a rough presentation of the 
simulation logic. It may be of some help to persons who attempt 
a detailed understanding of the program, but it is included here 
primarily to illustrate the gross program operation. 
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DEFINE FILE 10 (33 , 16, U, NEXT) ,1 (24, 4,U,NREC) ,2 (32 ,4 ,U,NREC) 
DIMENSION NOS (3, 4) , IP (3) 

WRITE (3,300) 

DO 1 J=l, 33 

READ (2,201)M, IQ, T, NOS 

WRITE (3,301) J,M, IQ, T, NOS 

1 WRITE (10* J)M,IQ,T,NOS 
MAXR = 24 

NUM = 1 
WRITE (3,303) 

11 DO 2 J=1,MAXR 

READ (2,202)KE, IP 
WRITE (3, 302) J,KE,IP 

2 WRITE (NUM* J)KE,IP 
NUM = NUM+1 

GO TO (99,3,99) , NUM 

3 MAXR = 32 
WRITE (3,303) 

GO TO 11 

99 CALL EXIT 

201 FORMAT (215, E5 . 0, 1215) 

202 FORMAT (415) 

300 FORMAT ('1 PART MASS CHG LIFE' 9X, 'DECAY PRODUCT PARTICLE NUMBERS 

1' / * 1 NUM TIME' 12X, * GROUPS OF THREE*/ ) 

301 FORMAT (' ' 14, 16, I5,E10.2, 4 (16,214) ) 

302 FORMAT (15, 19, 5X, 315) 

303 FORMAT ( ' 1 REAC THRESH RESULTING'/ 

1 ' NUM ENERGY PARTICLES') 

END 



Figure 2 
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The program PHYHE uses a short subroutine called BETWN which 
simply decides whether a number is below, within, or above some 
upper and lower limits • The listings for BETWN and PHYHE follow • 



SUBROUTINE BETWN (J,K,L,M) 

SUBROUTINE BETWN CHECKS WHETHER VALUE OF VARIABLE J LIES 
BETWEEN VALUES OF VARIABLES K AND L. 

M=1 

IF(J-K) 5,3,2 

2 IF(L-J) 4,3,3 

3 M=2 

GO TO 5 

4 M=3 

5 RETURN 
END 



Program : PHYHE 

DEFINE FILE 10 (33, 16, U, NEXT) ,1 (24, 4,U,NREC) ,2 (32 ,4,U,NREC) 
DIMENSION MO (3) ,IQ(3) ,MINM(3) ,MAXM(3) ,LP2(3,4) ,LP3(3,3,4) , 
IPN (15) ,NPR(3) ,TM(3) , JCS(3) 

C READ INITIAL EXPERIMENTAL CONDITIONS. 

1 READ (2, 201)M1,IQ1,KE, MO 
WRITE (3,301) Ml , IQl , KE , MO 
DO 3 IOUT=l, 3 

IF (MO ( IOUT) -137 ) 3,2,3 

2 WRITE (3,300) 

MO (IOUT) = 0 

3 CONTINUE 
NR=24 
IP=1 

IF (137+ (IQ1*M1) ) 12,5,12 
12 NR = 32 
IP = 2 

IF (496+ (IQl*Ml) ) 99,5,99 

5 MINE = KE - (KE/10) 

MAXE = KE 

DO 6 J=l, 3 

MINM(J) = MO( J) - (MO(J)/20) 

6 MAXM(J) = MO( J) + (MO(J)/20) 

INDX =1 

C FIND PARTICLES WHICH SATISFY SCAN CONDITIONS. 

DO 8 NRC=1, 33 

READ (10 'NRC)M,IC,T,LP2 

DO 8 J=l, 3 

CALL BETWN (M,MINM(J) ,MAXM(J) ,KEY) 

GO TO (8,7,8) ,KEY 

7 IPN (INDX) =NRC 
INDX= IN DX+ 1 

8 CONTINUE 
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C 



C 

C 



MAXNO=INDX-l 

FIND REACTIONS PRODUCING THESE PARTICLES 
KEY1=1 

DO 35 IB=1,NR 
KEY2=1 

READ (IP * IB) KET, NPR 

CALL BETWN (KET , MINE , MAXE , KEY) 

GOTO (35,21,36) , KEY 
21 DO 28 IL=i,MAXNO 
DO 28 IM=1, 3 

IF ( IPN (IL) — IABS (NPR (IM) ) ) 28 , 25 , 28 
25 IF (NPR (IM) ) 27,28,27 

BRANCH AT 27 WHEN PARTICLE NO. IM OF REACTION IB 
ONE OF THOSE SOUGHT BY THE SCAN SPECIFICATIONS. 

27 GO TO (41,42) ,KEY1 

28 CONTINUE 



IS IDENTIFIED AS 



C 



C 



WHEN NO PRIMARY PARTICLES SATISFY SCAN CONDITIONS CHECK TO 
FIND WHETHER DECAY PARTICLES ARE DETECTED BY SCAN 
KEY2 = 3 

DO 32 IM=1, 3 
KEY3 = 1 

NO = IABS (NPR (IM) ) 

JCS(IM) = NPR (IM) /NO 
IF (NO) 32,32,129 

129 READ (10* 10) M,IC,T, ( (LP3 (IM,MM,NN) ,MM=1, 3) ,NN=1,4) 

DO 32 IL=l,MAXNO 

DO 32 NN=1, 4 
DO 132 M2=l, 3 

IF (IABS (LP3 (IM,M2,NN) ) -IPN (IL) ) 132,30,132 
30 G T£ H (41 DECAY PR0DUCT LP3 SATISFIES SCAN CONDITIONS. 

130 KEY2 = 2 

GO TO (230,330), KEY3 
230 WRITE (3, 305) IM 
KEY3 = 2 



WRITE INFORMATION ON DECAY PARTICLES. 

330 KJ=IM 
I2=NN 
GO TO 146 
132 CONTINUE 
32 CONTINUE 

35 CONTINUE 

36 GO TO (97, 98) ,KEY1 

WRITE INFORMATION ON PRIMARY PARTICLES FOR SATISFACTORY REACTION 

41 WRITE (3,302) 

KEY1=2 

42 KJM = 0 

DO 44 KJ=1,3 
NO= IABS (NPR (KJ) ) 

IF (NO) 44,44,242 

242 READ (10 *NO)M, IC,T, ( (LP3 (KJ, II, 12) , 11=1, 3) , 12=1, 4) 

TM (KJ) = T 

JCS(KJ) = NPR (KJ) /NO 

JQ = IC*JCS (KJ) 

KJM = KJM + 1 
WRITE- (3,304) KJ,JQ,M 
44 CONTINUE 



0 



143 



146 



46 



47 



48 



49 



WRITE OUT DATA ON DECAYS 
GOTO (143,130,130) ,KEY2 
DO 50 KJ=1 ,KJM 
WRITE (3,305) KJ 
DO 50 12=1,4 
I1M=0 

DO 49 11=1,3 
ICS = JCS(KJ) 

MO (I1)=0 
IQ (Il)=0 

NUM=LP3 (KJ, 11,12) 

IF(NUM) 46,49,47 
ICS=NUM/IABS (NUM) * ICS 
NUM=IAB3 (NUM) 

READ (10* NUM) M, IC,T,LF2 
IF (M) 49,49,48 
MO (I1)=M 
IQ (II) =IC*ICS 
I1M = II 



OF EACH PRIMARY PARTICLE. 



CONTINUE 
IF (I1M) 232,232, 149 

149 WRITE (3, 306) TM (KJ) , (IQ (II) ,MO (II) , 11=1, IlM) 

232 GO TO (50, 132) ,KEY2 
50 CONTINUE 
GO TO 35 

97 WRITE (3,303) 

98 CALL EXIT 

99 WRITE (3,307) 

CALL EXIT 

201 FORMAT (615) „ ™™ m . , 

300 FORMAT ('0 ** PI MESONS ARE TOO NUMEROUS FOR A SCANNING REPORT / 

1 * THE SCANNERS ARE IGNORING YOUR REQUEST' ) 

301 FORMAT ('1 INITIAL EXPERIMENTAL CONDITIONS'/ 

' INCIDENT PARTICLE MASS='I5, ' MEV/C2 WITH CHARGE='I3/ 

vTHPmTP t'wi’Drv Qp PARTICLE BEAM = * 15 , * MEV * / 

■ i T*r 1 n iTn I tC 9 liUt 



1 

2 ' 0 KINETIC 
3 ' 0 SCAN IS 

102 FORMAT ('0 

10 3 FORMAT ('0 
1TIONS . ' ) 

t04 FORMAT (' 
t05 FORMAT ( ' 0 
106 FORMAT (' 
131X, 'CG = 

07 FORMAT ('0 



ENERGY — . 

FOR MASSES '15, ','15,', AND *15,’ MEV/C2 ' ) 

SUCCESS. RESULTS FOLLOW.'/) 

SCAN HAS DETECTED NO PARTICLES WHICH SATISFY YOUR CONDI 

PARTICLE ('II,') CG='I3,' MASS='I5) 

DECAY PRODUCTS OF PARTICLE ('ll,') FOLLOW.') 

DECAY TIME = ' E9 . 2 , ' CG = ' 12 , ' MASS = 

'12,' MASS = 'I5/31X, 'CG = '12,' MASS = '15) 

SORRY ' /' THE LABORATORY IS NOT ABLE TO PERFORM THE 



* 15/ 



1REQUESTED EXPERIMENT.') 
END 



COMPUTER SIMULATION OF A MASS SPECTROMETER 
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INTRODUCTION 



Every college physics department desires the best and widest 
range of equipment possible for its students to work with. For 
departments with reasonably large enrollments this ideal may be 
at least partially met. For departments with very small enroll- 
ments, however, the purchase of many pieces of major equipment 
cannot be economically justified. One answer to this problem is 
the simulation of physics instruments on a digital computer. The 
purpose of this paper is to describe how one such instrument, a 
mass spectrometer, has been simulated on the Anderson College IBM 
1620 computer, and how it is being utilized by the physics depart- 
ment. 



The program, MSSIM, is essentially a simulation of the Ealing 
Small Mass Spectrometer as described in Ealing Corporation^ 1969 
Teaching Catalog [also see J. W. Dewdney, Ameri can Journal of 
Physics 28 , 452 (I960)]. Ions from a source, with some thermal 
kinetic energy, are accelerated through a potential and deflected 
in a magnetic field to a detector which measures the ion current. 
The problem is to find the charge-to-mass ratio of the ions ob- 
served, and from this information the student tries to determine 
the isotope which is observed. 

The instructor inputs on punched cards the "sample" to be 
analyzed by the model. This sample consists of: 

1. The mass of each isotope present. 

2. The degree of ionization of each isotope. 

3. The relative abundance of each possible ion. 

He also inputs the value of the average initial kinetic energy of 
the ions, and initializes the random number generator to provide 
the "noise" current. 

The remaining variables are then input by the student on the 
console typewriter. These are: 

1. Strength of the magnetic field. 

2. Range of voltage to be swept. 

3. Increment to be used in sweeping the voltage. 

4. Width of the detector slit. 

The computer calculates an ion current for each value of vol- 
tage as it sweeps over the given range. The results are output on 
punched cards showing the particular voltage and the ion current 
associated with it. These cards may then be used as input to a 
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program which plots a graph of current versus voltage utilizing a 
CALCOMP plotter. 

Anderson College had no course in session during the semester 
in which this program was written for which it would have been ap- 
propriate to utilize this simulator, but this model will be used in 
the Modern Physics laboratory to be offered Semester I, 1969-70. 
Consequently, the Student Manual is not in its final form. In or— 
test the model under classroom conditions, however, an upper 
division student was engaged to perform two experiments with it. 

The purpose of the first experiment was to study dispersion 
and the results of using various detector slit widths. A plot of 
the expected graph of mass versus accelerating voltage was made 
and conpared to the experimental analysis of a sample containing 
93%K 39 and 7%K 41 , all isotopes singly ionized. By varying the 
slit width, different size detector current peaks were obtained 
(see Teacher’s Guide). 

In the second experiment the student was to identify an un- 
known element. An initial sweep was made to determine the loca- 
tion of the different ion peaks using a narrow slit width, and 
then localized sweeps were made of each peak using a wider slit 
width to determine the relative abundance of the isotopes present. 
The wider slit width was used in order to obtain flat-topped peaks 
which were more accurate to measure than the more pointed peaks. 

In this case the unknown was boron, and the sample consisted of 
18% singly-ionized B 10 , 1% doubly-ionized B 10 , 77% singly-ion- 
ized B* 1 , and 4% doubly-ionized B 11 . 

The results of these experiments proved to be quite satisfac- 
tory for both student and instructor. The first experiment led 
the student to a clearer picture of dispersion; the second gave 
him a better understanding of the use of the mass spectrometer as 
an analytical tool, as well as valuable experience in handling 
mass spectrometer data. Also, the student seemed to enjoy the 
experiment, especially the challenge of determining the unknown, 
which he did successfully. These conclusions are based both on 
the written report turned in by the student and conversations with 
him. 



The instructor was also satisfied with the outcome of the ex- 
periment and felt that the goal of giving the student some feel 
for using a mass spectrometer without a large monetary investment 
was achieved. This experiment showed that it was feasible to use 
this simulator within the normal three-hour period assigned to the 
Advanced Physics Laboratory course at Anderson College. In terms 
of actual computer time, these runs averaged four minutes per 
graph to calculate and 13 minutes per graph to plot. It is expec- 
ted that other experiments will be planned in the future to allow 
for even greater utilization of this simulator. 



STUDENT MANUAL 



Theory 

The mass spectrometer model to be simulated in our experi- 
ment is based upon the apparatus shown in Figure 1. The ions are 
produced in the ion source, and ejected with thermal kinetic en 
ergies ranging from 0 to 2AV electron volts. They are then 
accelerated through a potential of V volts to the source slits, 
which are assumed to be of infinitesimal width/ where they are 
collimated. They now have kinetic energy (neglecting thermal 
energies) 

y 2 rov 2 = qv 



where m,q,v are the mass, charge and speed of the ion. The ions 
are then deflected in a magnetic field, B. All ions that experi 
ence the same deflection and are brought together at an image 
point have the same momentum, given by: 



mv = aBr 



m 



where r = radius of curvature of the path followed by the ion. 
These two conditions are simultaneously satisfied by ions with the 
same charge-to-mass ratio. 



a = 



m 



2V 

B 2 r2 



thus the radius is given by: 

r = /2mV/B 2 q 



(3) 



(4) 



A radius of r = R 0 will bring the ions into the center of the 
exit slit, which has a width S. Ions passing through the exit 
slit will produce a measured current I d at the detector. 



The focusing of such an instrument is not perfect as there is 
an angular spread in the beam of any particular type of ion as l 
leaves the magnetic field, due to thermal kinetic energy. ^re- 
sult, the actual ion current reaching the detector is found by m 
tegrating over the width of the image slit, that is: 



I d = / J (r) dr + I n 

Slit 



(5) 



VJiAtVi 



where J(r) is the linear density dl s /dr of true ion current 
over the slit width as a function of radius of curvature, and I n 
is the random "noise" current of the system. 
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In this model, the component of source current due to each 
kind of ion present is determined by the following function: 



where Aj_ = relative amount of given ion in sample, Qi = relative 
charge of that ion, Iq = total amount of beam current leaving 
source slit, 5xl0~ 9 amperes. 

The ions will have a few eV of kinetic energy before they 
are accelerated. While this energy is best described by a one 
dimensional Maxwellian distribution, a triangular function is as- 
sumed in the model. It is further assumed that a dispersion of 
energy is equivalent to a dispersion of accelerating potential. 

By differentiating r in Eq. (4) with respect to V, it is pos- 

? 1 *1 I _ «•»_ • .1 1 _ 1 . 



The detector will measure the current caused by all ions 
falling between R 0 + S/2 and R 0 - S/2. This component of de- 
tector current is determined by the following: 



For S«.<R 0 then R 0 -S/2 < r < R 0 +S/2 gives radii of curvature 
for which ions will enter the slit. The integral in Eq. (8) is 
then evaluated for one of nine different cases, depending on which 




( 6 ) 



l AiQi 



i=l 



sible to find the difference in r, Ar , caused by a change, AV, 
in V : 




(7) 



Rn+S/2 
/ J(r)dr 
Ro-S/2 



( 8 ) 
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